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Abstract. We implement continued-fraction techniques to solve exactly 
quantum master equations for a spin with arbitrary S coupled to a (bosonic) 
thermal bath. The full spin density matrix is obtained, so that along with 
relaxation and thermoactivation, coherent dynamics is included (precession, 
tunnel, etc.). The method is applied to study isotropic spins and spins in a 
bistable anisotropy potential (superparamagnets). We present examples of static 
response, the dynamical susceptibility including the contribution of the different 
relaxation modes, and of spin resonance in transverse fields. 
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1. Introduction 

In the field of quantum dissipative systems one studies a subsystem consisting of a 
few relevant degrees of freedom coupled to the surrounding medium (bath), which 
has a large number of constituents (e.g., photons, phonons, electrons, nuclei) |lji2||3]- 
The (sub)system is not necessarily microscopic, but it can be a mesoscopic system (a 
Josephson junction, a magnetic molecular cluster, etc.) described by a few collective 
variables (phase across the junction, net spin) which under certain conditions can 
display quantum behaviour. Various fundamental problems can be addressed, like 
dissipation and quantum mechanics, decoherence, quantum Brownian motion, or the 
quantum-to-classical transition. The interaction with the bath, apart from producing 
dissipation, fluctuations, and decoherence, enables the system to interchange energy, 
momentum, and correlations with its environment and eventually relax to thermal 
equilibrium. For these reasons the study of open quantum systems is of interest in 
many areas of physics and chemistry. 

Classical open systems are handled as stochastic systems by means of Langevin 
and Fokker-Planck equations @]. This approach provides both a theoretical frame and 
computational tools, e.g., Langevin molecular-dynamics simulations. For few- variable 
systems, a powerful technique to solve Fokker-Planck equations is Risken's continued- 
fraction method 0] (a relative of Grad's moment approach for solving kinetic equations 
IS])- The non-equilibrium distribution W is expanded in a basis of functions and the 
coupled equations for the expansion coefficients Ct derived. An appropriate basis 
choice can render finite the coupling range (e.g., with the equation for Ci involving 
Ci-i, Ci, and C,;+i). Then these recurrences can be solved by iterating a simple 
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algorithm, which has a structure akin to a continued fraction. This method provides 
numerically exact results in problems of Brownian motion in external potentials, where 
closed form solutions are scarce (such classical problems are structurally similar to 
solve a Schrodinger-type equation). 

It results more delicate to deal with quantum dissipative systems First, 
phenomenological or non-standard quantisation poses problems with basic quantum- 
mechanical principles 0. Thus, one has to model the environment in a simple way 
(set of oscillators, 2-state systems, etc.), quantise it together with the system, and 
eventually trace over the bath variables. However, the resulting reduced descriptions 
are difficult to manage except in simple cases — free particle (or in a uniform field), 
harmonic oscillator, and few-state systems (e.g., S* = 1/2 spins). In general, the exact 
path-integral expressions for the (reduced) density operator g{t) are not easy to handle 
13 IH]; besides, the propagating function is highly oscillatory, rendering numerical 
methods unstable at long times j|9j. Quantum Langevin equations (Heisenberg 
equations for x and p including operator fluctuations) are of limited use beyond 
linear systems |10l 111) . Finally, under certain conditions (typically weak system-bath 
coupling), the density matrix obeys a quantum master equation [121 But again 
these equations can only be solved in a few problems. 

Due to their performance in classical systems (both translational E] and 
rotational ^HDj continued- fraction techniques were adapted to several problems of 
quantum Brownian motion in non-trivial potentials. This was done exploiting pseudo- 
probability representations of g |17| . Shibata and co-workers |18j applied continued 
fractions to solve master equations for isotropic spins {7i = —BzSz); Vogel and 
Risken tackled similarly quantum non- linear optical problems |19j . In phase-space 
problems, using the Wigner representation of the continued-fraction method for the 
Klein-Kramers equation was adapted [201 to quantum master equations of Caldeira- 
Leggett type [H] 12 1| (explicit recurrences were presented for polynomial and periodic 
potentials). As these approaches do not rely on the Hamiltonian eigenstructure they 
are applicable to demanding problems with (partially) continuous spectrum. 

Here we consider the following quantum-dissipative system: a spin in the magnetic 
anisotropy potential (H = —D Si — B ■ S) coupled to a boson (or bosonizable) 
thermal bath. For I? = we recover the familiar isotropic spin with its equispaced 
Zeeman spectrum; in a sense, the rotational equivalent of the harmonic oscillator. The 
anisotropy term —DS^ makes the problem tougher, say, the spin analogue of Brownian 
motion in double- well or periodic potentials. The environmental disturbances may 
indeed provoke a Brownian-type "reversal" of the spin, overcoming the potential 
barriers (Fig.P). In spite of the analogies, however, dissipative spin dynamics presents 
essential differences with translational problems, due to the underlying angular- 
momentum commutation relations [Si ,Sj] = ieijkSk |22l 12811^ . 

Our Hamiltonian describes paramagnets and superparamagnets — small solids or 
clusters with a sizable net spin {S ^ 10^-10''). For large S the physics is approximately 
classical (as in magnetic nanoparticles) and described by a rotational Fokker-Planck 
equation after Brown and Kubo-Hashitsume As the spin decreases, quantum 
effects come to the fore, as in magnetic molecular clusters where S ~ 5-25 |2Z]. The 
discreteness of the energy levels sensibly affects the thermoactivation processes, while 
the spin reversal may also occur by tunnelling when the field brings into resonance 
states at both sides of the barrier (Fig. An appealing dynamical description is 
given by master equations of Pauli type, for the diagonal elements of g ("balance" or 
"gain-loss" equations) JHEH!- These provide insight, while more refined treatments 
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Figure 1. Bistable energy levels of an anisotropic spin, Em = —Dm? — B^m, 
with S = 3 and _D = 1 at zero field (top left). An external field lifts the degeneracy 
£—m = Em (top right). Degeneracies are restored at multiples of D (lower panels). 



are less intuitive and difficult to apply. Nevertheless, to take into account coherent 
dynamics, like tunnel oscillations or the spin precession, one also needs off-diagonal 
density-matrix elements. But, as usual, solving the master equations for the full 
density matrix is not easy and several (often drastic) simplifications are required. 

In this work, following the spirit of Refs. TWTWW , we will solve master equations 
for non-interacting spins in contact with a dissipative bath by means of continued- 
fraction techniques. Exact methods available are limited to small spins {S < 5-10 
|29p. Our aim is to tackle arbitrary S, from the extreme quantum cases, S — 1/2 
and 1, to values as close as possible to the classical domain. This approach differs 
from those giving some continued-fraction expression for a certain quantity (typically 
relying on linear-response theory; see the review |30|1. in that here the full solution 
of the density-matrix equation is obtained by matrix continued- fraction methods. In 
this theoretical-computational frame we can study spins in a wide range of S {< 100- 
200) including their full dynamics: relaxation and thermoactivation, precession and 
coherence, as well as their possible interplay. 

The manuscript is organised as follows. We discuss the isolated spin and present 
the basic formalism in next section. In Sec. |31we introduce dissipative equations for 
a spin coupled to a bosonic bath, following the compact approach of Garanin et al. 
1^1221 with Heisenberg equations of motion for the Hubbard operators X™ — \n){m\. 
Master equations in the Markovian regime are discussed in Sec. 0] and fully specified 
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for several spin problems in Sec. \S[ In Sec. [HI we derive the chain of equations 
resulting from the perturbative treatment of probing fields (applicable to the non- linear 
response). In SecQwe manipulate the index-coupling structure of the density-matrix 
equations to obtain few-term recurrences suitable for implementing the continued- 
fraction algorithm. Numerous examples of its application to isotropic and anisotropic 
spins (superparamagnets) are given in Sees. |H1 and El we will check the results against 
exact formulae, whenever possible, and test heuristic expressions. We conclude with 
an assessment of our approach in Sec. 1101 and putting some auxiliary material and 
discussing technical issues in the appendices. 



2. Isolated spin (unitary dynamics) and Hubbard formalism 

Our starting point is a spin Hamiltonian j28ll33j including a magnetic anisotropy term 
and the Zeeman coupling to the field (units h ~ ks = gfJ-B = 1) 

H{S) = ~DSl-B-S . (1) 

This is the minimal model capturing the physics of (super)paramagnets. The term 
—D si has a bistable structure (for D > Q) with minima at Sz = and barrier top 
at ^2 = (Fig. Along with potential barriers (and degeneracies), an important 
consequence of the anisotropy is an unequally spaced energy spectrum ( [Appendix A| ).| 
Let us introduce the Hubbard (level-shift) operators 35, Ch. 1] 

X™s|n)(m|. (2) 

They form a complete system and any spin operator A can be expressed as 

A = EnmAnmX:^ , A„™ = {n\A\m) . (3) 

The expansion coefficients are simply the matrix elements of A in the basis defining 
the X^. (If not restricting ourselves to a multiplet with fixed S, we just need to 
introduce the corresponding indices \S m) and sum over them.) In particular, for 
the components of the spin operator one has Si = J2nmi''^\'^i\^) -^n' ■ Now, if we 
use the standard basis of eigenstates of and Sz, the required matrix elements are 
{n\Sz\m) — mSnm and {n\S±\m) — i^^n.m±i- Here S± = Sx i iSy are the ladder 
operators and = [S{S + 1) — m{m ± l)]^''^ the factors giving S±\m) — i^^\m ± 1). 
Then, the Si are represented by the single-sum forms 

m m 

In general, f{Sz) — X]m/("^)"'^m fo^' any operator function f{Sz); this gives the 
second-order "moments": S^ = and Sl + Sl = Y.JS{S -M) - m^] X-^- 

Concerning dynamics, the evolution in the Heisenberg representation of an 
operator is governed by i{dA/dt) — [A,H]. This plus Hamiltonian (QJ gives for X™ 

Ax™ = iA„„,X™ + 15+ {£+ - £- X™ i) + ii?- (Cn Xl^-' ~ £+ X™ i) (5) 

(see [Appendix C| ) where A„„j is the frequency associated to the m —> n transition 

A„,„ = e„ - e.m , £rn = - Dm^ - Bzm , (6) 

f Hamiltonian also describes a set of A'^ two-level systems interacting uniformly (Lipkin-Meshkov— 
Glick model) IM4I . The spectrum of 2^ eigenvalues splits in multiplets characterized by a certain 
S < N/2 and described by a pseudo-spin Hamiltonian Ti. = — D — BxSx', the excitation energy of 
each two-level element corresponds to the field and their coupling to the anisotropy parameter. This 
is a problem where the possibility of handling large values of S (large A'^) is important. 
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between levels of the diagonal part of the Hamiltonian. In the absence of transverse 
fields the evolution is simply = e'(*-*«)'^""'X™(to). In general, B± = B^± iBy 

couples the dynamics of the diagonal elements, X™, with the adjacent sub-diagonals, 
X^^^ and and so on. 

Finally, the density operator is expressed as g = J2nm QnmX™. The quantum- 
statistical average of X™ then follows from the trace formula {A) — Ti{gA) and reads 
(X™) — Qmn ( [Appendix C| |. This is important as it enables working with the density- 
matrix elements Qnm — ("-IpI'ti) or the X™ interchangeably, since all the equations we 
are going to handle are linear in X and their averaging thus trivial. 

3. Spin weakly coupled to a bosonic bath (dissipative dynamics) 

We address now the dynamics of the spin taking into account the coupling to its 
surroundings. Let us consider a total Hamiltonian including that of the spin 7i(S'), a 
bath of bosons Tih ~ X]q '^q'^ ^iid their interaction 

Htot - n{S) + Y.qCqFq{S) {a+ + a + ■ (7) 

Here Cq are coupling constants. The spin-dependent part of the interaction Fq{S) is 
typically a low-degree polynomial of S 28 , 36 . The coupling written is linear in the 
bath operators — only 1-boson processes are included, no Raman scattering involving 
two quanta (for interactions non-linear in the bath variables see Refs. [241 1311 IS7] '). 
Note that no counterterms are included in Titot to compensate for the coupling induced 
renormalization of the spin levels we will address this point below. 

The total spin-plus-bath system is unlikely to be in a pure state and a density- 
matrix description is required. For observables depending only on the spin, the 
required object is the reduced density operator g — Trb(£'tot), where one traces the 
bath out. For weak system-bath coupling a closed dynamical equation for g can be 
obtained by perturbation theory. This is the case of many problems in quantum 
optics, chemical physics, or magnetism ^ The equation has the generic form 
\{dg/At) — \H T g] -\- iR[g{T)\, where the relaxation term R adds to the Von Neuman 
evolution the effects of the bath. 

In the Hubbard framework, the Heisenberg time evolution of X™ = |n)(m| is 
governed by an analogous equation |^ 

dxr/dt = -i[xr,H] (8) 

The commutator generates the isolated-spin unitary evolution [Eq. ISJ] and i?™ 
accounts for the dissipation. When F{S) does not depend on the boson index [this is 
transferable to Cq; see Eq. (O], the relaxation term can be written as 

K-- f dr{/C(r - t) Fir) [F , X;"] - IC{t - r) [F , X,™] F(r)} . (9) 

This form is equivalent to the standard dissipative terms for g obtained by projection 
operators or cumulant expansions to second order O App. l.A]. The memory kernel is 
the autocorrelation /C(r) = {B{t-\-T)B{t))h of the bath operator i? — '^Zq^qi'^q^^ ~q)^ 
and reads K{t) = \cq\^[nq e+'"^-'^ + (n, + 1) e"'""^], with Uq = l/ie"^"/^ - 1) boson 
occupation numbers. This is how the temperature enters in the formalism, as the 
bath is assumed in equilibrium at the initial time r — )■ — oo. The operators without 
argument in Eq. Q are evaluated at t whereas F{t) — J2ki -^ki^Kr) introduces 
formally the previous history of the spin (cf. next section). 
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It is convenient to introduce the (coupling weighted) spectral density of bath 
modes J{uj) = X^qkqP tt6{uj — ujq). All quantities incorporating environmental effects 
can be expressed in terms of J{^)- For instance, the kernel /C(t) is given by 

JC{t) = / — J{lo) [n^ e+'^^ + {n^ + 1) e^'"^^] , (10) 
Jo 

with — l/(e'^/"^ — 1). A common functional form for the spectral density is 
J{ijj) (X Lo" (times a high-frequency cut-off at lod). The bath is called Ohmic when 
a = 1; this is realised by Kondo coupling to electron-hole pairs near the Fermi energy 
in solids (an example of bosonizable excitations from the ground state of a non-bosonic 
environment 6 ). For a > I the bath is called super-Ohmic; for instance, interaction 
with photons or phonons in three dimensions gives J{lo) cx 

We shall write J{uj) = Xuj" with A, determined by the |cqp, an overall measure of 
the coupling strength (classically, the damping parameter) . The characteristic "width" 
of the memory kernel, Tb, depends on the competition of 1/T and l/uJu, the bath 
bandwidth (the Debye frequency for phonons) . The relaxation term ||5J) was obtained 
treating the coupling perturbatively to second order for small Arb. 



4. Markovian (time-local) density-matrix equations 

Due to the integral term ^ the master equation is formally an integro-differential 
equation for X™. To second order in the coupling, however, the retarded dependences 
F{t) = J2ki ^ki^Kr) can be replaced by their unitary evolution, F{t) = U{T — t)F{t). 
Introducing these time dependencies in R only operators evaluated at t do remain. 

To illustrate, let us assume the Hamiltonian evolution simply given by X[(t) = 
gi(r-t)Afc, j^^^^^^ This can be plugged in i?™ and the resulting operator combinations 
expressed in the Hubbard basis and simplified using X^X^ — SmkXl^. 
This results in an equation of motion fully in terms of the X^{t), and linear in them. 
Carrying out these steps one actually gets the relaxation term ( [Appendix D.l| ) 

+{W:^„,+W^\^,)F„.r.F^mm' (H) 
— Snn' {J2e^e\m' FmiFim') ] A'™ . 

The coefficients include the matrix elements of the spin portion of the coupling, 
Fnm — {''^\F{S)\m), and the m n (complex) transition rates 

Wn\ra = W{^nm) . W^( A) = /jf dr e"'^'^ /C(t) , (12) 

evaluated at the level differences A„m = £n — £m- The form of the rate function T/F(A) 
emerges directly from the retarded dependences X{t) = e~'*^*~'^-'^A(i), which yield 
oscillating factors e~''^^ multiplying the kernel /C (r) in the integrand of Eq. On the 
other hand, the quantum-statistical average of the above i?™, using (A™) = Qmm gives 
an equation for the matrix elements {n\Q\m) which coincides what can be obtained 
from the standard, second order, relaxation terms for g [31 App. l.A]. 

Let us discuss when the conservative evolution of XI.{t) can be substituted 
by e'^'^"*^'^''' A^(i) in our problem. If one uses the basis of eigenstates of the full 
spin Hamiltonian (including the transverse terms |38l I39jl. such evolution is exact 
(then B± do not appear explicitly in the equation of motion, but only via A„m; 
[Appendix C| ). Further, if the transverse field is not too large, one can use the angular- 
momentum basis; then A™(r) ~ e'^'^~*-''^""'X™(i) gives the dominant Hamiltonian 
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dependences, providing an approximate relaxation term.§ This way of getting a time- 
local relaxation, without resorting to T — > oo approximations or semiclassical baths 
works when one explicitly knows the conservative evolution; apart from simple 
spin problems, it also applies to the harmonic oscillator |40t i41| . 



4-1. Relaxation term for couplings via S± 

In various important cases the coupling is realized through S± — Sx^iSy. For instance, 
F S± appears in Kondo coupling to electron- hole excitations and F ^ {Sz, S±} in 
magnetoelastic interaction with phonons [291 IS^ . Thus we will consider the form 

F{S)^V+{v{Sz),S-}+V-{v{Sz),S+} , (13) 

where {A, B} = AB + B A and rj± = r/^i: irjy are scalars incorporating the symmetry 
of the interaction (isotropic rj^ — rjy — 1; anisotropic rjx = \ and r^y = 0, etc.). The 
matrix elements Fnm = {n\F\m) are then | | Appendix D.2| | 

Fnm = Lrn,m-l5n,m-l + L*^+^^yn^n,m+l , Ljn.m' = ?7+['"(™) + vi'm')]£„i,m' , (14) 

where Lm,m' is an extended ladder factor with (im,m±i ~ [S{S+l)—m (mzbl)]-'^/^ ~ 
Although the operator v{Sz) commutes with 5*2, it modulates via [v{m) -\-v{m')] the 
matrix elements of S±, the ultimate responsible for transitions between the levels |m). 
The inclusion of a i^zS^ term in F does not lead to structurally new terms in the final 
equation and will not be considered here (it produces "dephasing" but not dissipation 
|31 Ch. 10]; in the language of magnetic resonance it modifies T2). 

The particularization of the relaxation term Hll|) to the coupling H13|) is done in 
[Appendix D.2| Invoking on it the secular approximation one is left with 

Rn ^ Ln,n-lL*„^.m-l (W^,t|n-1 + ^m|m-l) ^^^-l 
— ( |in,ri+lP ^n+l|n + |-^m,m+lP Wm+l\m 

+ |i„,„_ip I^„*_i|„ + W^„_i|„ ) (15) 

+ L*i^n+lLm,m+l (W'jtln+l + Wm\m+l) ^'^+1 ■ 

This corresponds to the rotating-wave approximation familiar in quantum optics, 
where counter- rotating, rapidly oscillating terms, are averaged out ( [Appendix D.3| ). 
Such manipulations seem not to pose problems for very weak coupling, while they 
simplify the treatment |42l I43j . Besides, the illustration of the continued- fraction 
approach will be cleaner disregarding non-secular terms in the master equation. 

To conclude, it is argued that the imaginary parts of the relaxation coefficients 
reflect a coupling-induced renormalisation of the levels, not genuine relaxation. In 
the bath-of-oscillators formalism this renormalisation is cancelled out by including 
suitable "counter-terms" in the starting system-plus-bath Hamiltonian "SJ. Here 
[ 311 I32| one cancels them by omiting the imaginary parts of Wn\mi redefining 
M^(A)^Re[/;"dre-'-^/C(r)]. || 

§ Classically, such corresponds to use R = —Alls' x (S x zB^^) as relaxation term in the 
Landau-Lifshitz equation, fully keeping the Hamiltonian precession dS/dt = {S X B^g) + R. (Here 
B^ff = —&H./dS\ the form S X &H./dS follows from dSi/dt = {Si,W} via the Poisson brakets 
{Si ,Sj} = eijkSk 1341 '). The 2-component of B^g retained in R incorporates the anisotropy field 
Sa ~ 2DS, the dominant energy scale in superparamagnets. 

II The formalism takes as previously assessed whether such a renormalization is physically meaningful 
for a given coupling, and if so (e.g., the Lamb shift), it is considered to be already included in Ti.(S). 
This has the advantage of making of Ti. the experimentally accesible Hamiltonian, instead of the bare 
one which may be difficult to determine. 
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4-2. Elements and structure of the resulting density-matrix equation 

The basic ingredients of the master equation will be the energy differences An^, 
transverse fields B±, ladder factors, couphng matrix elements {n\F\m) ~ im,m', and 
transition rates Wn\m = W{Anm)- AH properties of the bath enter via the rate 
function, which can be expressed in terms of the spectral density J(w) ( [Appendix E| ) 

W{A>0)^J{A)nA, W^(A < 0) J(|A|)(n|A| + 1) . (16) 

As — 1 / {e'^/^ — 1) boson absorption and emission rates are related by the detailed 
balance condition W{A) = e"'^/^ 1¥(-A). This ensures, under certain conditions, 
convergence to the Gibbs distribution at long times |42l HH] . 

The rates in the relaxation term (|15ll involve adjacent levels only, Wm\m±i = 
W^(Am mj-i), while i?™ connects X™ with AT^ differing in indices by at most 1. Thus 
in the sequel we will compactly write 

which together with Eqs. (O and (jHl, gives the working equation 

It is worth mentioning that the density-matrix equation was obtained within a, 
though approximate, fully quantum treatment, not introducing any phenomenological 
relaxation or assuming preconceived structures for the equation. 

Note finally that if the transverse field is set to zero B± = 0, the diagonal part 
of Eq. (|18|l becomes a closed system of balance equations for the level "populations" 
Nm = AT™, as in the Pauli master-equation approach 

Nm = + Wm + W+ N^+l . (19) 

As Amm = 0, the Hamiltonian part does not show up and the dynamics is purely 
relaxational (W,- = >V";,"Ii, >V,„ = W:^,;^, W+ = W";;^+i). For example, in the 
case of an isotropic spin in a static field we can always choose the z axis so that 
Ti. = —BzSz- However, if we want to study resonance phenomena, we need the full 
equation (|18|l to account for transverse probing fields. Further, even if B± = (or 
when using exact eigenstates of the full Hamiltonian) , the simple balance structure (|19|l 
is broken by terms like X^^l or AT™^^ if not resorting to the secular approximation 
( [Appendix D.2| ). For these reasons we will focus on the full density-matrix equation. 

5. Density-matrix equations for specific spin-bath problems 

Here we particularize equation (jl8|l to an isotropic spin with truly linear coupling 
F ~ S± , as that to electron- hole excitations (an Ohmic bath) , and to anisotropic spins 
with F ~ {Sz, S±}, corresponding to interaction with phonons (a bath from Ohmic to 
super-Ohmic depending on the lattice dimensionality). In the classical limit, the first 
coupling yields field-type fluctuations in the spin Langevin equations j^SJI^SI, whereas 
the second produces anisotropy-type fluctuations O EH) (the spin analogue of 
force- type and frequency- type fluctuations in mechanical systems j46j). 
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5.1. Density-matrix equation for isotropic spins 

Let us consider a spin ?i = B S, with the linear coupling F = i(r;+S'_+?7_5+). This 
corresponds to a constant v{Sz) = 1/4 in Eq. (|13ll . Then L.„i,m' = \'n+^m,m' in the 
matrix elements F„m and hence Ln^n'L*m,m' = j\'n\'^^n,n'£m,m' , with jyyp = ?7+f7-- On 
the other hand, as the transition rates involve adjacent levels Wj^^j^^i = W{Am^rn±i)i 
and here Am,m±i = ±-6^, only two rates appear. Introducing wq = Wra\m-i for the 
m — 1 ^ m transition (decay for Bz > 0), and W^m-i|m = wqC'^^^^ by detailed 
balance, the relaxation term (|15|l reduces to 

K = u;o{C-i4„-i^,ri^-i [(^^+4)+e~n^Li+4-i)]^"+e-^44„^r+Y} (20) 

Here we have introduced y = B^/T, the single factor £,n — [S{S + 1) — m{m + 1)]^^^ 
(then ^+ — £m and £^ ~ tm-i), and assumed isotropic coupling rj^ ^ rjy = 1 
(|?7p simply rescales wq). Plugging this relaxation term and the Zeeman transition 
frequency A„m = —Bz(n — m) into the master equation l|18() . one finally gets 

X™ = -is,(n-m)x;" + + 

+ z«o{C-iC-iXrY - 1 [(^n + 4) + e-^(^Li + 4-i)]^r + e-y£nimX::^+'} (21) 

This equation corresponds to the master equation of Garanin |31j (he also included 
2-boson Raman processes). For B± = the diagonal elements obey the balance 
equations ^ with W„ = wo^,2„_i, W,„ = -■u;o(€^ + e-2'£^_i), and >V+ = wne-^ll^. 

In an Ohmic bath J{uj) — Xuj the rate function is W{/S) = XA/{c^/'^ — 1) (this 
form is valid for both signs of A; see [Appendix E| |. Then wq = XTy/{l — c^^) which 
in the limit y — *■ goes over the classical rotational-diffusion constant wq ^ XT 
(independent of -Bz); thus A from J(w) plays the role of the Landau-Lifshitz damping 
parameter. Actually, taking the limit S* — > oo in the balance equations 03 App. A] 
one gets the classical Fokker-Planck equation (in a longitudinal field), and the 
correspondence is established as 2 AS' All- 

5.2. Density-matrix equation for superparamagnets 

Next we consider anisotropic spins, H = —DSz—B-S, with a coupling linear in S± 
but with Sz-dependent "coefficients", as it occurs in spin-lattice interactions. There 
F = ^{Sz, r?+S_-hr?-S'+}, corresponding tow (Sz) = 52/2 in Eq. lO. For jryp = 2 this 

gives Ln,n'L*m.m' = \^n,n' £m.m' ■, with the modulatcd factors 4n,m' = (m + m')£m,rn'- 

Then the density-matrix equation H18() goes over 

+ \£n-l£m-l{Wn\n-l + W.m\rn-l)X^_i -\- \in(-m{Wn\n+l + Wm\jn+l)X^n+l 

- i(^lW„+i|„ +^lW^™+i|,„ +^l_iW^„-i|„ +^l„iW^™_i|™)X™ , (22) 

where A„m = —[D{n -t- m) -\- Bz\{n — m) and we have introduced the corresponding 1- 
index notation = (2m -f- l)£m ( [Appendix B| ). This equation was derived in Ref. [32 
for the study of the archetypal magnetic molecular cluster Mni2. The replacement 
(■m ~* (not affecting the Hamiltonian part) accounts for the 5z-dependent coupling 
and results in an extra m dependence of the relaxation term with respect to the case 
F S±. It can be seen as a level-dependent "damping" Xcfi{m) ^ X{2m -|- 1)^, 
decreasing as the anisotropy barrier to ~ is approached fFig. IBljl : a spin analogue 
of position-dependent damping in translational Brownian motion. 
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The rates Wn\m involve adjacent levels but no simplification arises here due to 
the non-equispaced spectrum of anisotropic spins A„i.m±i = ±[D{2m ± 1) + B^]- 
To get the rate function W^(A) one can compute the distribution of bath excitations 
J{ijj) = '^q\cq\'^ T^6{uj — bjq) with a Dcbyc phonon model Wfcs = \k\ and replace 

— !■ 'Y^^ Jd'^k, integrating over wave-vectors and summing over branches. For 
magneto-elastic couphng one has |cqp ~ ujq so that d'^fc x |cfcsp ~ x |fc| gives 

spectral densities J oc to" evolving from Ohmic a = 1, for phonons in one dimension, 
to super-Ohmic J (x for d = 3. The corresponding relaxation functions are 

W{A) = A A"/(e'^/^ - 1) , VA . (23) 

This unified form, instead of Eq. (|16|l . is discussed in [Appendix E| Note that phonon 
velocities Vg, coupling constants, etc., are subsumed in A from J(w) = Aw". 

6. Response to probing fields: perturbative density-matrix equations 

With the master equations one can describe the non-equilibrium evolution from one 
stationary state to another. A system can be made to "relax" either by subjecting 
it to a "force" (a magnetic, electric, stress field, etc.) or by removing it after having 
kept it for a long time. Then the question is how the infusion or withdrawal of energy 
is shared by the system's degrees of freedom. Alternatively one can apply a force 
oscillating with frequency fl; this provides a time scale 1/ft whose competition with 
the intrinsic scales of the system permits to analyse its different dynamical modes 0H| ■ 

To reflect intrinsic properties the probe should be suitably small, not altering 
the nature of the studied system. This has the advantage of allowing the use of 
perturbation theory in the treatment. In this section we will replace B hy B -\- 6B in 
the spin density- matrix equation, and treat it perturbatively in SB, getting a chain of 
coupled equations. Each level will be tackled sequentially with the continued-fraction 
treatment of SecQ giving the spin response to the perturbation (susceptibilities). 

To alleviate the notation, we first write the density-matrix equation (|18|l in the 
following compact form (including all n and m) 

X = iA{B)X + W{B)X . (24) 

Here iA X stands for all the Hamiltonian part (we put the "i" to remind us of this) 
and W X for the relaxation term. The field enters via iA linearly [Eq. ^] and non- 
linearly through Am,m±i = i[D{2m±l)+Bz] in the rates W{A) ~ A" ua- Let us now 
augment B by a probing field, B B -\- b{t) u, with u = [ux, Uy, Uz) a unit vector 
along 5B and h{t) — bcos{nt) its magnitude. To obtain linear susceptibilities, we 
expand field-dependent parameters / — f{Bx, By, B^) to first order: / ~ /o + /i h{t), 
with /k = (d''//d6'')/K! derivatives with respect to the amplitude d/d6 = 'YiUidsi- 
For the coefficients of the master equation this gives iA — iAo + \Aib{t) and 
yV ~ Wo + yVib{t) (the former is exact as A is linear in B). 

Although modulated quantities like iA and W have the parametric time 
dependence f{t) = f[B -f h{t)] ( [Appendix D.3| ), our dynamical variable X{t) does 
not need to evolve as some function of i? + h{t). Thus, we seek for a solution of the 
form [no (t) in b] 

X{t)~Xoit) + Xi{t)b. (25) 

We compute now iA x X and W x X to first order, replace them in the dynamical 
equation (|24l) . and equate terms with the same power of b, getting (X_i = 0) 

X^ = (iAo + Wo)X^ + (iAi + Wi) cos(m) . (26) 
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The perturbative structure is clear: original equation with unperturbed coefficients 
(first term) plus their derivatives ( • )i times the previous order result. Thus the lower 
level X^_i acts as a forcing (source) term when solving the equation for X^. 

To get the long-time stationary response^ when all transients have died out, we 
introduce Fourier expansions (subindex for order in b, superindex for the harmonic) 

X.{t) = ^Y^.X^:^^"'''' . (27) 

and go order by order. The harmonics e"^'^* generated at each n will coincide with 
those of the forcing cos(rit) X^_i(t), because Eq. 1)26(1 is linear in and itself 
is not multiplied by oscillating terms (additive sources). The zeroth order equation 
Xq = (iAo + Wo)Xo has no sources. Then, only v = Q is left in the Fourier series, 
Xq = -X"q°\ and the static response obeys = (iAo + Wo)-X'o°''. At first order the 
forcing is cos(ili) Xq — i(e+'^* + e"'^*)^!,"'' and only the harmonics v = ±1 get 

excited {X^°^ = 0). Thus we plug Xi = i(X^^^e+'"* +xj"^^e-'^*) into Eq. ^ and 
equate Fourier coefficients at both sides, getting the remaining equations 

/ = (iAo+Wo)x(°) + 

\ ±il7Xr^ = (iAo+>Vo)xl±^) + (iAi+Wi)x(°) ^'^^ 

These equations are to be solved sequentially by the continued-fraction method, with 
the previous order acting as a forcing on the next. 

Some final remarks. Technically, we use an equation- of -motion approach to 
obtain the response, not relying on Kubo-type linear-response-theory expressions '48|. 
This allows proceeding systematically to higher orders to get non-linear suceptibilities 
(harmonics of the excitation generated by non-linearities). Then, one finds terms of the 
form cos'(rit)XK-i in the equation for X„, due to the non-linearity of the 

relaxation term [which includes T1^(A) = AA"/(e^/^ — 1)]. In quantum Brownian 
motion as described by the Caldeira-Leggett equation |^, due to the high-T 
approximations plus Ohmic bath (corresponding to W ~ AT), the relaxation term 
does not depend on the system potential: R = —i'y[^{x—y){dx — dy)-\-T{x—y)'^]Q{x,y) 
OT R — ■y dp{p-\-T dp)W{x,p) in the Wigner representation. In particular, R does not 
depend on the forcing, and such higher-derivative terms vanish {Wi>2 = 0). They 
are also absent in classical spins and dipoles |49[ EO] , where the relaxation term does 
depend on the field, R = —All S x {S x Bcs), but linearly. Anyway, as we have seen, 
the yV;>2 terms do not affect the calculation of the linear responses. 



7. Continued-fraction methods for spin density-matrix equations 

To solve the master equations we will cast them into the form of 3-term recurrence 
relations suitable to apply the continued- fraction method 0]. This is related with 
schemes of solution by tri-diagonalization, like the Lanczos algorithm or the recursion 
method |51| . In Brownian motion problems it shares elements with the expansion into 
complete sets (Grad's) approach for solving kinetic equations jS]. The non-equilibrium 
distribution W is expanded in a basis of functions (Hermite or Laguerre polynomials, 
plane waves, spherical harmonics, etc.) and the partial-differential kinetic equation 
is transformed into a set of coupled equations for the expansion coefficients Ci. 
Approximate solutions can be obtained by truncating the hierarchy of equations at 
various levels. But as a matter of fact, to obtain manageable equations the truncation 
needs to be performed at a low level. In the continued-fraction variant, instead of 
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truncating directly, one seeks for bases in which the range of index couphng is short 
(ideaUy, the equation for Ci involves Cj_i, Q, and C^+i). Then these recurrences 
between the C'i can be solved exactly by iterating a simple algorithm, which has the 
structure of a continued-fraction.^ 

This technique was exploited to solve classical Fokker-Planck equations for few- 
variable systems in external potentials jl4lll5l[T?>j . Compared with direct simulations, 
continued-fraction methods have several shortcomings: (i) the basis choice is quite 
problem specific, (ii) the stability and convergence of the algorithm may fail in some 
ranges of parameters, and (iii) it does not return trajectories (which always provide 
helpful insight). When the method works, however, its advantages are valuable: (i) 
no statistical errors, (ii) special aptness to get stationary solutions (long times), (iii) 
high efficiency, allowing to explore parameter ranges out of reach of simulation, and 
(iv) the obtaining of the distribution W (also insightful) . 

This, together with the lack of quantum Langevin simulations, motivated several 
adaptations of the continued- fraction approach to quantum problems |18[ 1191 I2(jj . 
The master equation was transformed into Fokker-Planck-like form using pseudo- 
probability representations W of the density matrix, such W expanded in appropriate 
bases, and the recurrences for the coefficients derived from the kinetic equation and 
solved by continued fractions. As this approach is not based on the Hamiltonian 
eigenstructure, it is invaluable for systems including continuous parts in the spectrum 
(e.g., Morse and periodic potentials). Notwithstanding this, for systems with discrete 
levels only, the density-matrix equation for gnm already has an index-recurrence 
structure and such transformation-expansion protocol (often cumbersome) may be 
bypassed. This is our purpose in this section (cf. Refs. |5HI 1541 Sec. V] for S — 1/2). 



7.1. Index- coupling structure and vector 3-term recurrences 

Let us begin writing the master equation fll^ compactly as X™ = J2n'm' Q™'n^ -^n/ ' 
with n' = n — 1, n, 7i + 1 and to' = m — 1, m, m + 1, and with coefficients 

^m,m— 1 >*,m.m— 1 ^m,ra i r> p— ^m,m-t-l 

Qm,m— 1 i. R f~ r^ni^m \ \ i ")/\n7i,m y^m,m + l J. R 

n,n ~ 2 m ^n,n ~ ^ '-^nm C K„ „ Wn,n ~ 2 + m \'^^ I 

n,n+l = ^ '^n,n+l ^ ^2 n '^n,n+l — ^^n,n+l 

The matrix associated to this linear system has dimensions {2S + 1)^ x (25* + 1)^. For 
small spins it can be solved directly. However, already at = 3 the size is 49 x 49, 
becoming 441x441 for S" = 10 (Mni2 or Fcs), and rising to 1156x1156 for S = 33/2 
(Feig) and to 2704x2704 for S = 51/2 (Mnas). Thus, if one is tempted to study 
mesoscopic spins in this way, let alone pursue the classical limit, soon faces large 
matrices. 

The problem gets simpler if B± = 0, as the system splits into uncoupled 

recurrences inside each sub-diagonal X'^+^ = These 
can be solved separately by scalar continued fractions, as in the approach of 
Shibata |18| . Nevertheless, in problems involving coherent dynamics the diagonals 
are typically coupled and such a strategy is not applicable. They remain coupled even 
when B± = (or using the exact eigenstates of the full Hamiltonian) if not resorting 
to the rotating- wave approximation [Eq. ljD.3|l ] . 

1 For a nice brief survey of the relation between ordinary series expansions, orthogonal polynomials, 
recursions, and continued fractions, see Ref. 1521 . 
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This calls for more generic methods. Our aim is to retain the spirit of Shibata's 
approach by converting the 2-index recurrence X™ = J2 Q^n" , into a 1-index 
vector recurrence. To this end, we first rewrite the equations highlighting the index- 
coupling structure 





\rm — l 


/nm,m— 1 
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x^m — 1 










x^ - 
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Forgeting for a moment about the upper indices m', we see that the equation 
for Xn^ only involves X^li (first column), Xn^ itself (second), and X^]_i (last). 
Thus, if we build up "vectors" c„ including all X™ for a given n we will have 
a 3-term recurrence between them, with some matrix coefficients Qn,n'- That is, 
Cn = Qn,n-ic„_i + Qn,nCn + Q„_„+i c„-|_i . This could be tacklcd by matrix continued 
fractions, replacing the original (25* + 1)^ x (25* + l)'^ problem by one with 25* + 1 
steps but with matrices (25 + 1) x (25 + 1) ( [Appendix G| ). This approach actually 
converts the ^ S'^ problem into an overall ^ S"^ one and reduces significantly the 
storage demands. Compared with previous exact techniques, this will allow us to 
increase the possible values of S dramatically (up to ^ 100-200). 

Explicitly, the required (25"+ l)-tuples c„ and (25'+ 1) x (25+ 1) matrices Qn,n' 
comprise (with averages in (c„)„j = {X™) one deals directly with gnm) 
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and 



Then, introducing the custom notation Q„ = Qn,n-i, Qn = 'Q 
the density-matrix equation H3()(l is transformed into the 3-term canonical form: 

n = -5,---,5. (32) 



It is worth emphasizing the simple attainment of the sought recursion form in 
comparison with the Fokker-Planck-like approach with pseudo-distributions (as a 
consequence of the discrete spectrum). On the other hand, if the original equation had 
included couplings with X™^^ would have arrived at 5-term vector recursions (e.g., 
for biaxial spins, or not invoking secular approximations). Nevertheless, recurrences 
involving more than 3 terms can be "folded" into 3-term form by introducing 
appropriate block vectors and matrices | | Appendix D.3| ) . Alternatively, if the elements 
breaking the 3-term structure are suitably small, they can be treated iteratively 
(avoiding to enlarge dimensions), in a way similar to the forcing terms below. 
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7. 2. The forcing terms ( iterative calculations ) 

The equations kiVLX^ ~ (iAg + Wo)^k + (iAi + >Vi)X„„i of Sec. El can be 
manipulated analogously, as they inherit the structure of the original master equation. 
The main novelty are the forcing terms (iAi + Wi)X ^-i, which will be handled here. 

Recall that bold letters in Eqs. (|27|I - H28|I standed for all elements, to which 
we will have to add indices for perturbative order and harmonic. Then to avoid too 
baroque expressions we introduce the simplified notation Y = x''^^ and Z = X^'-p . 
The 0th order equation (|27|l can then be written as [cf. Eq. H18|l ] 

+ (Wo)™r-7' + (Wo)™;r + (Wo):i:r+t' , (33) 

where the index stands for absence of probing field. The first line corresponds to 
IAoXq'^'' and the second to WqX'^^\ Similarly, the first order Eq. (|28|l reads 

ifi Z„ =iA^j^Z„ + ^ijjj. ^ ~l^^Z^_^ + ^B^_{(.^Z^ —i'^Z^j^-^) 

+ [(Wo)":"7' z^^-i + im)Z:n K + {y^^)n;n+V ^r+Y] + .C(^) • (34) 

Again the custom terms stand for iAp-X'^^'' and Wo-X^i^' while the sources /™ account 
for (iAi + Wi)x''q\ This is determined by the previous order result Y — Xq"' 
and the field derivatives of iA and W. Using 6B ~ b {ux,Uy,Uz) and A„m = 
— [D{n + m) + Bz]{n — m) we obtain 

+ [(Wi)"„"i7' y:^,' + {w,r:;r y:^ + Y:r+V] , (35) 

where u± = ± iuy and (Wi)™;^"' = d(>V™;T')/dS^, because W only depends on 
Bz (recall that our relaxation term is approximate in B± ; Sec. 0)) . 

Now, to convert the (25 + 1)^ x (25" + 1)^ systems and (|S3l into vector 
recurrences, we proceed just as in Sec. 17.11 for the parts KiilX^ ^ (IAq + Wq)Xi^, 
while we introduce appropriate forcing vectors /„. This gives 

Q-Cn-l + QnCn + Q+c„+i = -/„ , Q„ = Q„ - ^117 I . (36) 

Now {cnU = or and (/„),„ = /™, while (QnUrw = Q™^' and (QtUm' = 
Q^'^i as above. The modified central matrix Q„ (I is the identity) incorporates the 
left-hand sides KiflX,^. The source terms (absent for k — 0) can also be written as 
/„ = dbQ,^ c^'r/'' + dbQn ci"'~^'' + db(}+ c':^+i \ with c^""^'' the previous order resuh 
and the matrices differentiated with respect to the probing field: ( • )i = d( • )/d6. 

Equation H3t)|l has the canonical form permitting to apply directly the continued- 
fraction algorithm of [Appendix G| Besides, apart from perturbatively, the form (|36|l 
also follows from the original master equation c„ — Q^c„_i -I- Q„c„ -I- Q+c„+i 
through Laplace transformation (for t-independent Qn)- Then one just identifies 
Qn = Qn — sl (i.e., Kifl s, the Laplace variable) and /„ = Cn{t = 0) ^ p(0). This 
would allow tackling initial-value problems (not forgetting the cautionary remarks of 
[Appendix D.3[ ). 

7.3. Spin response and susceptibilities 

Once the recursions are solved we have all density-matrix elements (c„)m = (X^) = 
Qmn and any observable can be obtained from the trace formula (A) — J^nm^nmAmn- 
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For instance (Si) — X)nm(^l'^*l'^)('^")"i' which connects directly the spin response 
with the continued-fraction results. To get explicitly the response to SB = 6cos(f2i) 
we insert the expansion X™ (t) ~ + 1 (^™e+ + Z^e~ ) into the above average 

[cf. Eqs. and (123, Y and x[^'^ -> Z] 

(S^) = Enmi^\S^HY^ + | { ("I 1"^) e+ + [E„™ ("I 1^) Z„"] C' } . 

Comparison with (Si) — {Si)o + ^(xe+'^* + x* e~'^*) gives the static response (S'i)o 
and the dynamical susceptibility x(^)- terms of its real part x' and imaginary 
part x" the time-dependent response is {Si){t) — {Si)o = b{x'cosi}t + x"sinrit). 
Proceeding to higher orders, as sketched in Sec. El Xg"-* x[^^ X^^^ xf^ ^ • • ■, 
the non-linear susceptibilites (higher harmonics) would follow similarly: x'^'^H^) ~ 

8. Application to isotropic spins 

Now we proceed to apply the discussed methods to solve the density-matrix equations 
for various systems. We will start with Garanin's master equation for isotropic 
spins, a. = B ■ S, with a simple linear coupling to the bath F — t] ■ S. The bath is 
assumed Ohmic, J(w) = Aw, with rate function W{A) = AA/(e^/^ - 1) (Sec. EHJ. 
There are several analytical results for the static and dynamical response of isotropic 
spins, which will be used as benchmark for the continued-fraction approach. 

8.1. Matrix coefficients of the recurrences 

Comparing the relaxation term (|2()|l with the generic form (|17|l we identify the 
relaxation coefRcients W^^T' of the isotropic spin. Plugging them in Eq. (|29|l for 
the coefficients Q„ and the Zeeman level differences A„m = — (n — m)Bz, we have 

{^m.m— 1 n n / ^m,m — 1 p. 

Qn.n-1 = tn-lim-l \ Qn'n+1 = 

Q™;™! = -(i/2)s+c-i (Q)+ QZ':r+i = -m)B.e„ (3?) 

r Q™r"' = (i/2)i?-c-i 

Qn<Q™r - -i{n-m)B,-{wo/2)[i£l+il)+e-yiel_,+£l_,)] (38) 

Recall that y — B^/T, the decay rate is wq — XTy/{l — e^^), while im = (-ti and 
= €~ with £^ = [S{S + 1) — m{m ± 1)]"^^^. From these coefficients one also gets 
the derivatives dbQ„ for the treatment of probing fields, B ^ B + h{t) u, completing 
the specification of the 3-term recurrences H32|l and (|36|l . 

8.2. Thermal- equilibrium response 

Analytical results. The statics of isotropic spins Ti. — —B^Sz can be studied in full 
analytically, giving us the opportunity to test the continued-fraction solution of the 
master equation. The magnetisation = (Sz) is given by the Brillouin function Bs- 

mz = SBs{0, Bs{0 = {^ + ^)cth[{l + ^)(]-^cth{^i) , (39) 
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Figure 2. Magnetisation = (Sz) vs. field of isotropic spins. Left panel: 
ruz vs. Bz/T for various small S. Right panel: Reduced magnetisation m^/S 
vs. f = SBz/T for increasingly large S. The symbols are continued-fraction 
calculations and the lines Brillouin functions 1391 [the dashed line is the classical 
limit L(^) = cth(^) - 1/5]. 
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Figure 3. Temperature dependence of the susceptibilities of an isotropic 5 = 10 
spin in various external fields. Left: reciprocal longitudinal response 1/X||; the 
lines are Curie— Brillouin susceptibilities 1401 . Right: transverse response; lines 
Xi_ = rriz/Bz [Eq. I41i ] . The symbols are continued-fraction results. 



with the scaled field variable C = SB^/T Sy). For SB^ ^ T, Bs I and 
saturation ~ S" is reached. The longitudinal susceptibility x\\ = d{Sz) / dB^ follows 
using (cthx)' = 1 — cth^a: and Bb^ — dy/T: 

XII =5(5+l)/r-{(5+i)2cth2[(l + 5L)^] -i cth2(5L^)}/r. (40) 

Expanding the hyperbolic cotangent as cthx ~ l/x + x/Z, the second term goes over 
-2S{S + l)/3r as Bz 0, completing the Curie law xo = S{S + l)/3r for the 
initial (zero-bias) susceptibility of isotropic spins. On the other hand, the response to 
a probing field perperdicular to Bz is given by the transverse susceptibility: 

XX = mz/Bz , rUz = {Sz) . (41) 

Here we used Van Vleck's formula to get x± = Xxx with the transverse-field 
dependent energy levels xu = Sm['^(^«^'")^ ~ e"^^™ where diS = de/dBi. 

For SBz <C T, one has niz — BzX S{S -\- 1)/3T, recovering Curie's law from this side. 
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Numerical results. Figure [21 shows the agreement of the continued-fraction resuhs 
with the analytical {3^)- The curves for small spins exhibit saturation to the 
corresponding S at large fields. Increasing S, up to S* = 100, we can follow the 
evolution towards the classical Langevin magnetization / 5* — > cth (^) — 1 [Eq. (|39|l 
with the leading terms 1 + 1/25* ~ 1 and cth(^/2S') ~ 25/^]. At 5 = 20 the result 
is already close to the classical asymptote. However, this depends on the field range 
observed, for quantum behaviour is found whenever the discreteness of the energy 
levels is important, and this can be attained by increasing sufficiently Bz- 

The agreement of the magnetisations, which is numerically exact, ensures 
agreement for the susceptibilities. Still we have computed X|| — dmz/dBz 
directly from the equilibrium fluctuations of the spin X|| = {{Sz) — {Sz)^)/T, to 
check the proper continued-fraction obtainment of second-order moments (SiSj) = 
^^^(n|5'iS'j|m)(c„)m. Figure 01 shows l/x|| vs. temperature for a moderate spin. 
In a small Bz there is a straight-line dependence in almost all the range [Curie law 
Xo^ = 3T/S{S + 1)]. Raising the field we observe deviations upwards (maximum in 
X||) at sufficiently low T. This occurs because x\\ is the slope of the magnetisation 
curve and, at high Bz/T, saturation S takes place and the slope drops to zero. 

Figure also displays the transverse susceptibility. For a quantum spin x± is not 
easily expressed in terms of averages in the absence of perturbation, so we resorted to 
applying directly a small transverse field and computing x ~ {Sx) / B^- The agreement 
with Eq. H41|l is remarkable (recall that the relaxation term is approximate in the 
transverse field. Sec. 0J. We see how x± is reduced as T increases, approaching the 
Curie regime. At low T the magnetization saturates, ~ S*, and the curves tend to 
the constant values XJ-(^ = G) = S/Bz- 



8.3. Dynamical response 

Let us turn now to the dynamics of isotropic spins. Wc will consider the response to 
probing fields h cosiVlt) parallel and perpendicular to the bias field Bz . 



8.3.1. Analytical results. For 5B \\ B, on replacing Bz ^ Bz + fecos(r2i) in 
the balance equations (|19|l (with the coefhcients of Sec. 15.111 and expanding to first 
order in b one obtains equations determining the longitudinal susceptibility. They 
form a system of few coupled equations for small spins and can be solved analytically 
I2I1I13EH1. For = 1/2 one finds the Debye form [(•)' = d( ■)/dy; y = Bz/T] 

X||(f^)-^JT^, r = wo{l + e-y), (42) 

where ruz = ^th{^y) and m'^/T is the equilibrium response. The characteristic 
relaxation time is t = l/F. As the decay rate is wq = W{—y), with W{y) = 
XTy/{ey - 1), one has F = W{-y) + W{y) (see | AppenduTEl ) . Next, for 5 = 1 
the susceptibility can be written as PT| 

<^ViA2+jf]r^ _ Fi 

^11^ ^ T (Ai+ir!)(A2+ir!) ~ TFi+if!" ^ ' 

+ For the statics we use a weak spin-bath coupling A = 10~^ in the density-matrix equation I21i . 
We know that in the absence of transverse field its diagonals are decoupled (after the rotating wave 
approximation). Then detailed balance iy(A) = e~^/^ W(— A) ensures that the Gibbs distribution 
is its statiorary solution. Thus the static continued-fraction results must be independent of A. 
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Figure 4. Dynamical susceptibility X\\(^) for ^ = 1/2 (left panels) and 
5=1 (right). The temperature is T = 1, A = 10~^, and results in various 
fields are shown. The lines are the analytical expressions 1421 and 1431 and the 
symbols continued-fraction results. The different panels are the real part x' (top)i 
imaginary part x" (middle) and Cole-Cole plot x" vs. x' (bottom panels). 

Here Ti = r(2chy -f- l)/{chy + 2), with T = W{~y) + W{y) again, and Ai,2 are 
the non-zero eigenvalues of the matrix associated to the system of balance equations 
Ai_2 = r[2ch(ii/) ^ l]/ch(i?/) (Aq = corresponds to the thermal equilibrium 
solution). This formula can be expressed as the sum of two Debye terms (cf. Eq. (|48|l 
below). But as Ai is numerically close to Fi, the susceptibility is nearly single Debye 
with an effective relaxation time t ~ l/Fi. 
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Figure 5. Imaginary part of the susceptibility for S = 1/2, 1, 2, 3, 5, 20. 
Again T = 1 wliereas tlie fields are chosen to fix ^ = SBz/T = 0, 2, 4, 6 for 
the different S. The solid lines are the single-relaxation time approximation 1441 
and the symbols continued-fraction calculations [normalised by the Curie xo = 
S{S + 1)/3T]. For 5 = 6 we also display numerical results for S = 50 and 100 
(rhombi), indicating convergence to a classical limit. 



For arbitrary S the longitudinal relaxation comprises 25-1-1 modes "SH" (the 
rank of the matrix associated to the balance equations). In general their amplitudes 
tti and rates = Aq < Ai < • • • < A2S need to be obtained numerically j47|. 
This makes it difficult to derive closed expressions for the susceptibility, x{^) = 
X^i=i'^i/{^ + if^/Ai), but still something can be said on the relaxation time. An 
overall measure is provided by the integral relaxation time, ri„t = /q dtSM{t)/SM{0), 
where 5M{t) is the response to a small field change SB^ at i = 4 . Its advantage is 
that one can by-pass the computation of the time constants and amplitudes |44II45[IT7) 
by getting Tint directly from the low-frequency susceptibility xi^) — x(l ~ i^^Tint + ■ • .) 
which can be obtained in closed form (see |Appendix F] for explicit expressions for Tint). 

Finally, to have some analytical formula for the susceptibility at arbitrary S, one 
introduces a single-relaxation time approximation 

•^«'"> = ^TTik;^ 

This is possibly the most popular expression in the modelling of relaxation 
experiments. By construction this heuristic form is correct at low frequencies. The 
accuracy for arbitrary f2 will have to be assessed in each problem addressed. 
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8.3.2. Longitudinal response. Figure ^ displays susceptibility spectra for a spin 
5 = 1/2 in a number of fields, showing the agreement between the analytical and 
continued-fraction results (again numerically exact). The curves go down as the field 
is increased. This is mainly due to the reduction of the equilibrium part of the response 
m'^/T [Eq. H42|l ]. which is the slope of the magnetisation and decreases with increasing 
Bz. The peaks of the imaginary part x": E^nd the maximum slope of the real part 
x', occur at n ~ r. Then, as the relaxation rate F = W{—y) + W{y) increases 
with Bz ( [Appendix E| |, the curves shift to higher frequencies. Finally, as the response 
comprises a single Debye factor x ~ 1/(1 + i^:^)) plotting imaginary vs. real parts 
(Cole-Cole or Argand representation) gives perfect semicircles jST,. 

The longitudinal response for S" = 1 is also shown in Fig. 0] In comparison with 
spin one-half, the susceptibilities are higher and more sensitive to the bias field. This 
can be related with the magnetisation curves of Fig. |21 which for 5=1 have larger 
slopes and change faster with Bz (the parameter coupling to the field, i.e., 5, is larger 
now). On the other hand, as the deviation of Eq. H43|) from a single Debye is small, 
the Cole- Cole plots are nearly semicircular. 

For arbitrary 5*, finally, we compare the continued-fraction results with the 
heuristic formula = X||/(l + i^^^'int)- Figure El shows that the agreement is 

very good in general, implying that the relaxation modes A^ are quite grouped (on a 
logarithmic scale) . There are only small deviations at the peaks in intermediate fields 
(^ ^ 2-4), in accordance with Garanin's findings in Ref. On the other hand, 

this figure (and Fig. ISJ illustrates that in order to reach the classical asymptotes, the 
quantities need to be scaled appropriately.* 

8.3.3. Paramagnetic resonance. We conclude with the response to transverse 
probing fields. These incorporate S± in the Hamiltonian, not commuting with the 
dominant ^z-dependent part, and provoking transitions \m) — > |rn ± 1) between the 
unperturbed levels. These transitions result in peaks at the frequencies /S.m,m±i = 
£m — Emiti in the imaginary part of the susceptibility @HI (absorption line-shape). 
Classically the phenomenon corresponds to the matching of the oscillating field with 
the Larmor precession of the spin. Finally, as the quantum transverse response 
{S±) = '^„i(-m{^m±i) involves ofi^-diagonal elements of g (coherences), one refers 
to the precession in transverse fields as coherent dynamics. 

As the level spacings of an isotropic spin are all equal Am,m±i — ^Bz, the 
phenomenology for the different S is qualitatively similar. Figure El shows the 
susceptibility x_l = Xxx for S = 1/2, for which we have exact formulae 



[To get this equation one has to solve the full density-matrix equation (|21|1 : for 
S = 1/2 the form coincides with that obtained from phcnomenological Bloch equations 
identifying T2 — l/r2.] The continued-fraction results and Eq. H45() agree to all 
computed figures. In particular, the numerical results duly fulfill the basic relation 
x(— f^) = [x(ri)]*, yielding even x'(^) ^nd odd x"{^)- The imaginary part shows 

* We used the scalings — > SB^/T, — > m^/S, and x — * x/^i^ + !)■ In general, we increase 
S keeping the amount of Zeeman energy and the anisotropy barriers constant (and hence finite at 
5 — ► 00). For the Hamiltonian H = —DSl — B-S this amou nts to fix D and SB, while 
introducing more levels with S (their spacing decreases as A ~ 1/5, [Appendix X} . Correspondingly 
I58II47I App. A], frequency and damping are scaled as H oc 1/S and A oc 1/5 (recall that A ^ All/25). 
Alternatively one can use the scaling combination Q/A (or some Qt), as we did here. 
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Figure 6. Transverse susceptibility x±{^) of isotropic spin S = 1/2 in a 
magnetic field Bz = 0.1. The lines are Eq. 1451 and the symbols continued- 
fraction results. Left: real and imaginary parts for A = 0.01 at T = 0.5. Right: 
detail of the effects of the damping (halved to A = 0.005) and of the temperature 
(doubled to T = 1) on the absorption line-shape x'j_{^)- 



peaks at \Q\ — B^, the level separation, accompanied by zig-zag with sign change 
of the real part. Decreasing the spin-bath coupling A the absorption peaks become 
narrower and higher, as in a forced and damped oscillator. Here, this can also be 
attained by changing T (right panel). These behaviours, captured by the line- width 
T2 in Eq. (|45|l . reflect the "smearing out" of the energy levels due to the bath coupling. 

With this simple example we have introduced the basic phenomenology of 
magnetic resonance and some factors influencing it. On the other hand, the perfect 
agreement of the continued-fraction results with the exact solution lends confidence 
in our handling of the non-diagonal elements of g, required to compute {S±) = 
12m^m9m,m±i- This is important for the subsequent application to spins in the 
anisotropy potential, where there are less analytical expressions to compare with. 



9. Application to anisotropic spins (superparamagnets) 

Now we consider spins with a Hamiltonian Ti. — ~D — B ■ S and a quadratic 
coupling F ~ SzS±, motivated by the spin-lattice interaction in (super)paramagnets 
1221 EHl- Correspondingly, we use the rate function VK(A) = AA'^/(e^/^ — 1) of a 
a = 3 super-Ohmic phonon bath. At variance with the uniform Zeeman spectrum, 
the anisotropy results in non-equispaced levels (Fig. ^ and hence in several rates 
Wm,\m±i = W{Am,m±i)- Thus this problem is a spin analogue of translational 
Brownian motion in non-harmonic potentials. 

9.1. Elements of the density-matrix recurrences 

To solve the density-matrix equation (|22() by continued fractions we convert it, as 
explained in Sec. 17.21 into a vector 3-term recurrence of the form Q~c„_i -I- Q„c„ -t- 
Q^c„-|_i = — /„, with {cn)m = Qmn- The matrix coefficients Q„ comprise Hamiltonian 
and relaxational contributions, which for this problem follow comparing Eqs. (|22|l 
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and H3U|) [we include A„m = — [D{n + m) + Bz]{n — m)]: 
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Recall that the replacement im ^ = (2m+ in the relaxation parts comes from 
5*2 in '--^ {Sz, S±}, leading to "position-dependent" damping. The field derivatives 
of these coefficients give the source terms /„ for the treatment of probing fields 
(see [Appendix E| for dW/dBz); then (Q„ Q„ — k'iUI. With the density-matrix 
recurrence (|36f) so specified we can now apply the continued-fraction algorithm of 
[Appendix G[ to solve it. 

9.2. Thermal- equilibrium response 

Again we begin discussing briefly the static properties in a longitudinal field. Compact 
expressions for the magnetisation to^ = m e"'''^™ follow for small spins. 

Introducing d = D/T and y = B^/T, we have —Psm = dm? -\- ym. Then for spin 
one-half — ^(e''/^ — e~^/^)/(e^/^-|-e~^/^), equal to the isotropic result (as D enters 
in the two levels equivalently) . However, for S — 1 and 2 the magnetisations read 

2e'^shy _ 2e'^shy + 4e4'^sh(2y) 

^ l + 2e'ichy ^ l + 2e'i chj/ + 2e4'^ ch(2y) ' ^ ' 

which are valid for both D > (easy-axes anisotropy) and D < (easy plane; then 
the energy levels of Fig. ^ are turned upside-down). Notice that the to = level does 
not contribute to X^m™^"^"^™ ^lut contributes "phase space" in Z = X^m^"^'^"- 

For D > the states with m — zLS have the lowest energies in weak fields and 
are only separated by 2BzS. Then the magnetisation curves have the convex features 
of the isotropic-spin case (Fig. [T] open symbols). In contrast, for D < the curves 
depart from zero slowly (exponentially); the low- field ground state is then to = 0, well 
separated from the first excited level (by \D\—Bz). Indeed, for 5* = 2 and D < 0, when 
the field makes m = 1 the new ground state, the magnetisation is again stabilised at 
niz — 1 until it "jumps" to rUz — 2. The jumps become steeper as T is decreased. 

The longitudinal susceptibility for S" = 1 follows differentiating niz in Eq. 146() : 
the transverse response x± can be obtained from Van Vleck's formula fSec. l8.2|l . The 
results for the initial (zero-bias) susceptibilities are 

12 2 l-e~^/^ 

^" ^ T 2 + e-^/T ' ^^^:D 2 + e-W^ ' = " ^^^^ 

At high T both expressions recover the isotropic susceptibihty: X|||d=o — X±\d=o ~ 
2/3T [here 3(3 + I) = 2]. At low temperature one has x\\\d^oc> = l/T for D > 
(2-state like) in accordance with rUz ^ thy ~ TOz|s=i/2- Thus in both limit ranges 
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Figure 7. Magnetisation curves of anisotropic spins with S = I (left) and 5 = 2 
(right) both for positive and negative D atT = 0.2. Lines, Gibbsian formulae I46i : 



symbols, continued-fraction results (with A : 
3 



2.5 
2 

1.5 
1 

0.5 




10-3; vd. footnote in Sec.lOl. 







X||(D 
X_i(D 
S(^+1 


= 1) . 

=-1) = 

= 1) o _ 

=-1) . 

m — 




v.. 












S=1 
6^=0 



















0.5 



1.5 T 2 



Figure 8. Equilibrium susceptibilities vs. temperature of 5 = 1 spins with D = 1 
(circles) and D = —I (squares). We show longitudinal responses (filled symbols) 
and tranverse (open), together with Eqs. 1471 (thin lines). The thick solid line 
is Curie law, approached by all curves at high temperatures, and the dashed line 
the Ising asymptote [for X||(^ > 0)]- 



X|| obeys 1/T laws, with the factor 1/(1 + ^e^^/-^) governing the intermediate T 
crossover between the isotropic-spin and Ising regimes. In contrast, for D < Q, the 
longitudinal susceptibility goes to zero exponentially at low T. Again this is due to the 
m — ground state for easy-plane anisotropy. Finally, as for the transverse response 
at low T, it tends to the constant limits x^U^oo = ^/D and X-lU^-oo = 2/lZ?!. 

Textbook examples of the longitudinal and transverse susceptibilites |55| are 
displayed in Fig.|Sl showing the agreement between the analytical expressions and the 
continued-fraction results (the longitudinal obtained as X|| = {{Sz) ~ {Sz)'^)/T and 
the transverse from x±(^) using a small il/X = 10""^; cf. Sec. 18.2(1 . This agreement, 
together with the magnetisation curves of Fig. [7| indicates that we are handling 
properly a quantum system with non-equispaced levels, as well as its transverse 
response, which requires off-diagonal density-matrix elements. 
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9.3. Dynamical response 

9.3.1. Analytical results. Again exact expressions can be obtained for small 
S from the balance equations ((T^ . with coefficients W" — Wm\m-i^m-iJ — 
-(^mW^m+i|m + and W+ = For g onc-half the coupling 

model considered does not produce relaxation (see [Appendix Bj l. The first non-trivial 
case is S" = 1, whose longitudinal susceptibility comprises two Debye factors 0^ 

Here A^ are eigenvalues of the balance-equations matrix 

Ao = o, Ai,2 - (r+ + r_) T V(r+ - r_)2 + aw+w^ , (49) 

with the rates r± = 2[Vl^(A±i_o) + A±i^o)] and transition probabilities w± = 
P±i|o = 2 W{A±ifi) [recah that W{A) = A A°/(e^/^ - 1), A±i,o = -{D±B,), and 
i±ifl = iv^]- The amplitude a G [0, 1] controlling the weights of the two summands 
involves Ai^2 and Aofr = (w+ -I- w^)/Zm'^^ the initial decay rate of the magnetization 
|T7] . As Zm'^ = 2(chy-|~2e'')/(2ch2/-j-e~'^) we have AoH' ^ Ti as D/T ^ 0, recovering 
the susceptibility of isotropic 5 = 1 spins [Eq. (|43|l ] . We are not aware of exact results 
for X|| (ri) of larger spins or for the transverse dynamical response with D 0. 



9.3.2. Longitudinal response. As the range of parameters that can be explored 
is wide, we concentrate on low temperatures (and eventually, weak fields), the 
experimentally most interesting range in superparamagnets. It is convenient to 
introduce reduced anisotropy and field parameters 

a^\D\S''/T, ^^SB,/T, h^^/2a. (50) 

The latter is Bz in units of 2|£)|S' which is of the order of the anisotropy field at the 
minima or the field for barrier disappearance | |Appendix A| ). As mentioned before, 
when comparing different S we scale parameters keeping a and ^ fixed. 

Let us begin with S = 1. Its dynamical susceptibility is shown in Fig. Inland the 
full agreement between the analytical and numerical results. The curves evidence two 
relaxation modes [(25* + 1) — 1, the equilibrium Ao = 0] . Inspecting the structure of the 
corresponding eigenvectors [47| . the low- frequency mode, Ai, can be associated to over- 
barrier crossings and the faster mode, A2, to transitions between neighbouring levels 
(±1 0, intra- well dynamics). The over-barrier process dominates the response at 
weak fields; the intra-well is active but by symmetry its contribution to (Sz) practically 
cancels out (but not to (5^), the Kerr relaxation observable). Increasing Bz the 
spectrum losses the potential barrier at Be = D {he = 1/2) and the fast transitions 
between adjacent states take over. For Bz>D the two modes are still separated, 
because the levels are not equispaced yet. Finally at high enough fields a Zeeman 
spectrum is approached and the isotropic susceptibility recovered [Eq. (|43|l ]: then the 
two modes are close (in frequency) and x(^) approaches again a single Debye form. 

Let us now address the response of larger spins. Figure ITUl shows that, although 
we should be finding 25* modes in x(^)j they appear gathered in two main groups: 
the over-barrier mode Ai and a bunch of high-frequency modes, related to intra-well 
transitions. This, in turn, leads to a phenomenology akin to that of 5 = 1. In this 
figure we have chosen the Argand plot (x" vs. x') where competing modes are resolved 
in two neat semicircles. They evolve into one in the limits of low and high field. In the 
low Bz regime the response is dominated by the over-barrier dynamics and there is a 
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Figure 9. Longitudinal susceptibility spectra of an anisotropic S = I spin with 
D = I and A = 10"^ at T = 0.1 and various h = Bz/2DS, below and above 
barrier disappearance Be = D {a = 10, f = /i • 2(T = — > 18). Lines, exact 
two-mode Eq. I48i : symbols, continued-fraction results; the arrows trace the field 
evolution. The susceptibility is normalised by the equilibrium H — > value. 
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Figure 10. Cole-Cole plot of the longitudinal susceptibility of spins 5 = 3,6, 12. 
The anisotropy parameter D is adjusted to fix cr = D /T = 10 and the fields 
such that i = SB^/T = 0, 2, 4, 6 (/i = 0, 0.1, 0.2, 0.3). The lines are single-Debye 
approximations and the symbols continued-fraction results. The susceptibilities 
are scaled by xo = + 1)/3T and the y axis range is half the x axis one. 



good agreement with a single Debye form x{^) = x/(l + with r ~ 1/Ai. But in 
contrast to = 1, the intra- well modes are clearly manifested at fields h* well below 
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Figure 11. Dynamical susceptibility x\\ vs. temperature of an anisotropic spin 
5" = 10 with D = 0.5 at several frequencies (in units of A/D^). Real parts (solid 
symbols), imaginary parts (open); single-Debye Eq. I44i (lines). The numbers by 
the SI = 10~* curves correspond to the regimes discussed in the text. Considering 
D and T given in Kelvin, the values used are close to those of Mni2. 

the field of barrier disapearance he — 1 — 1/25' [Eq. (|A.5|l ]. Indeed, in the classical limit 
it was estimated that above only h* ~ 0.17 the fast modes can significantly compite 
with the over-barrier process j59| due to the thermal depopulation of the upper well 
(cf. Fig.nj. For large S (say S'> 10) we obtain h* approaching such classical result. 
Besides, we see that the onset of the intra- well modes depends on the spin value. That 
is, h* = h*{S), increasing for decreasing S. This seems natural because the results 
should recover h* ~ 0.5 as S ^ 1. Equivalently, at a fixed h the semicircle of the fast 
modes (the left one) is less developed the smaller the S is (panel h = 0.2). Eventually, 
at large a single-Debye again describes x(^) for a-U S. (For further discussion on 
h*{S), the modes interplay, and analytical approximations, see Ref. |47|.') 

9.3.3. Superparamagnetic blocking. Finally we present the results at low fields, 
the range most studied in nanomagnets (2DS ~ lOT in Mni2) in a way closer 
to experiment [501 There one varies T at a fixed il, because the exponential 
dependence of the relaxation time r cx exp{AU/T) permits to span various decades 
in rtr in an easier way. The results show the phenomenology of superparamagnetic 
blocking — a maximum in the magnitude of the response at some intermediate fi- 
dependent temperature (Fig. 111(1 . This is different in nature from the maxima 
exhibited by the equilibrium x\\ for Z? < 0, or for > in an external field (Figs. 13 
and 13). The dynamical blocking is due to the competition of r with l/Jl and the two- 
fold rolex played by T. It unblocks the over-barrier transitions at low temperatures, 
enabling the spins to follow the oscillating field bcos{ilt), but sufficiently high T also 
provokes the thermal misalignment of the spins from the field direction, degrading the 
response. 

Let us follow the process in some detail. (1) At low temperatures, r ^ l/fi, 
the probability of over-barrier crossings is negligible and the dynamics consists of 
transitions at the bottom of the wells (with a small averaged projection onto the 
field). (2) Increasing T the spins appreciably depart from the minima, and the response 
starts to rise with T. However, as the thermo-activation is not efficient enough, the 
response {Sz){t) ^ 6 (x' cos fit -I- x" sin fit) sizable lags behind the field, as manifested 
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Figure 12. Continued-fraction results (lines) for x±{^) of an anisotropic 5 = 10 
spin with D = 0.5 and X/D^ = 3 ■ 10~*. The calculations were done at zero field 
and T = 10 (cr = 5). Real part, left panel; imaginary part, right. Vertical lines: 
loci of the transition frequencies Am^m+i = D{2m + 1), m = 0, 1, ■ ■ ■ , 9. 



by the large x" ■ (3) At higher temperatures the over-barrier mechanism becomes 
more efficient; the reponse continues increasing, but becoming more in-phase with 
the excitation (x' dominates). (4) If T is further increased, however, the thermal 
agitation also provokes misalignment of the spins from the field direction. Then the 
magnitude of the response exhibits a maximum and starts to decrease; this occurs at 
a temperature Tb such that T{Tb) ~ l/fi. (5) Eventually, at high T, the spins quickly 
adjust to the equilibrium distribution corresponding to the instantaneous field. Then 
x' goes over the equilibrium susceptibility (ex 1/T) while x" drops to zero. Note finally 
that, in agreement with the previous subsection, this low- field phenomenology is well 
described by a single Debye form, as Fig. 1111 shows. tl 

9.3.4- Paramagnetic resonance of anisotropic spins. We conclude with the spin- 
resonance behaviour of quantum superparamagnets Recall that a field oscillating 
perpendicular to the anisotropy axis provokes transitions |m) ^ |to ± 1) between the 
unperturbed levels. Computing the response along such field {S±) = '^rn^rnQin,m±i 
requires off-diagonal density-matrix elements and falls outside a Pauline master 
equation for the populations gmm- 

The induced transitions result in peaks in the absorption line-shape x"(^) 
at the frequencies Am,m±i. For an isotropic spin all level differences were equal 
Am,m+i = Bz. Here, however, the anisotropy yields non-equispaced levels Am,m+i = 
D(2m+l) + Bz and one would expect multiple peaks in x"(r2), with the corresponding 
zig-zags in the real part x'{^)- At zero field the 25-1-1 levels are degenerated by pairs, 
m with —m (Fig.^, and we should find only S peaks at the locations = D(2to-|- 1). 
The largest frequencies (A^ ~ 2DS) correspond to transitions at the wells (|r7i| ^ S), 
while those near the barrier top appear at low O (|m| ~ 0, Ab ^ D). Finally, as we 
saw in the isotropic S' = 1/2 spin (Fig.O, the absorption peaks have finite width and 
height due to the damping and the temperature. 

fj For low Q we have found numerical instabilities of the continued-fraction results with very large 
AU/T. Some accuracy problems also arised in isotropic spins at very large B^/T. They can be 
attributed to the exponential dependences of the relaxation rates, giving tiny numbers at certain 
critical places; this is known to happen in the classical case under the same limit conditions. 
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Figure 13. Effects of the external field, damping, and temperature on the 
line-shape x" (^) °f anisotropic 5 = 10 spins with D = 0.5. Left: results for 
A/D^ = 3 ■ 10"* and cr = 5 at = (as in Fig.IEJ, = 0.25 (= D/2) and 
= 1 {= 2D; cf. Fig. nj. Right: reference zero-field line of left together with 
results halving the damping at the same T, and halving T while keeping A. 

Figure El shows these features for S — 10. Starting from the rightmost x" 
peak, associated to the ground-state transitions, the intensity of the peaks decreases 
with 17, as they progresively involve transitions between higher levels, thermally less 
populated. On the other hand, the peak width is not uniform in f2. This is due 
to the S'2-dependent spin-phonon interaction, F ~ SzS±, which gives an extra m 
dependence of the relaxation term (compared to the coupling F ^ S±). This enters 
via the modified ladder factors — {2m+l)'^£'^ and can be seen as an effective level- 
dependent "damping", Xce{m) ^ A(2m+1)^. Therefore the transitions between upper 
levels (lower m and f2) correspond to a reduced effective damping, and those peaks 
become narrower and higher. Overall, the competition of the thermal depopulation 
and the reduced damping yields peak heights initially decreasing as fl is reduced and 
rising again at low frequencies. ff 

We conclude with various effects on the line shape x"(il). The application of 
a field lifts the degeneracy by pairs and the peaks split (Fig. left). At the first 
resonance ~ D the levels become degenerated again (m = and m = — 1, m = 1 
and 771 — —2,. . . ; Fig. ^ and the energy differences, A = 0, 2D, AD,. . . , are just 
half-way those of zero field, A = ID, 3D,. . . , D{2S — 1). Then the "side peak" to the 
right merges with the one to the left of the neighbouring peak (curve not shown) . For 
Bz > D the peaks split again (they simply crossed) and at the first even resonance 
Bz = 2D they merge again but on the original locations (the spacings correspond to 
the Bz = ones plus A = D{2S-l)-\-2D from the ground state; Fig. HI). FigureHHalso 
shows the sharpening of the peaks when decreasing the damping (as in the isotropic 
S = 1/2 spin, with the new feature of non-uniform widths) and the dramatic reduction 
with temperature of the intensity of the lower Q lines. Those transitions involve higher 
levels whose thermal population gets exponentially reduced as T is lowered. 

ft Note that W {Am,m±i) can also add to the dependence on m )58l . except for an Ohmic bath at 
high T (then W ~ AT). However, the bare factor = S{S + 1) — m,(m -\- 1) does not. It is 
geometrical, giving the factor (1 — z^) in the classical Fokker-Planck equation |47| which accounts 
for the increased phase space at large angles z = cos 1?. 
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10. Summary and discussion 

In the field of quantum dissipative systems one usually works with the reduced density 
operator g of the subsystem of interest (tracing over the bath). In several problems, 
due to the weak system-bath coupling, one can derive perturbatively a closed equation 
of motion for g — a quantum master equation. This plays the role of the Fokker- 
Planck kinetic equation in classical problems. We have addressed one such system, a 
quantum spin with arbitrary S* in a dissipative environment, and solved exactly the 
corresponding density-matrix equation by implementing continued-fraction methods. 
As the density matrix is obtained in full, coherent dynamics is included along with 
the relaxation and thermoactivation. 

The continued-fraction method belongs to a family of exact methods in condensed 
matter and statistical physics and has been fruitfully exploited in problems of 
Brownian motion in external potentials. Here wc took advantage of the index- 
recurrence structure of the density-matrix equation to bring it in the form of a 
few-term recurrence relation, suitable to apply continued fractions. For simple spin 
problems this had been done by Shibata and co-workers, exploiting the decoupling 
into independent equations for the diagonals, gm,m+k ~ F{Qm',m'+k)^ and solving 
them by scalar continued fractions. In general such decoupling does not take place 
(e.g., in the presence of transverse fields), and matrix continued fractions are required. 
This has been the contribution of this work, allowing the obtainment of numerically 
exact solutions of the full density-matrix equation for quite generic spin problems. 
Besides, compared with previous exact techniques, the spin values affordable have 
been increased significantly (up to 5* ~ 100-200 on an old laptop). This range of S 
should be enough for studying the evolution to the classical limit in many problems 
(one of the central issues in open quantiim systems). Reaching large 5* is also important 
when effective spin models arc used to describe interacting 2-level systems. 

Technically, we have worked within a Hubbard formalism (Heisenberg equations 
of motion for the operator basis X™ = |n)(m|), whose main advantage is being 
compact. This is not essential, however, and was used as intermediate step; the 
equations for X™ are linear and their averaging gives directly the standard equations 
for the density- matrix elements Qnm = {n\o\m). We have focused on stationary 
responses, for the obtainment of which this numerical method is specially suited (in 
contrast with path-integral propagation schemes affected by sign problems at long 
times). Our starting point was to convert the time-dependent master equation into a 
perturbative chain of stationary density-matrix equations with each step solvable by 
continued fractions. We have worked it in full for the linear dynamical susceptibility, 
but the extension to get non-linear responses is systematic. Besides, upon Laplace 
transformation, a number of time propagation problems could also be tackled (of 
the type "evolution between stationary states", not "system-meets-bath" problems 
requiring time-dependent coefRcients) . 

From the outset the implementation was done in its general form (with matrix 
continued fractions), not taking advantage of the splitting into diagonals of simple spin 
problems. This made the initial tests tougher while it allowed proceeding smoothly to 
more general problems, not enjoying such decoupling. The implementation has been 
simpler than in a Fokker-Planck-type approach with pseudo-distributions, as it avoids 
the transformation of the density operator into some phase-space representation, 
the expansion in complete sets of functions, and eventually the manipulation of the 
coefficient recurrences (as in quantum Brownian motion in phase space). Similarly, the 
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implementation has been easier than in the continued-fraction solution of rotational 

Fokkcr-Planck equations for classical spins and dipolcs, as some aspects are simplified 
in the quantum case. For instance, the finite number of discrete levels results in 
finite continued fractions, and convergence or termination problems are fortunately 
by-passed. Thus, some numerical instabilities found are to be attributed to accuracy 
problems when handling tiny numbers; actually they appeared in parameter ranges 
already problematic in the classical case (very low T and fl). On the other hand, the 
finite number of steps in the algorithm can be carried out by hand for small spins. 
The approach can actually be called semi-analytic, which is the reason behind the 
numerically exact agreements found with explicit solutions. 

In this frame and with these tools wc addressed the statics and dynamics of 
spins with arbitrary S in contact with a thermal bath. We have considered the 
familiar isotropic spin, H = —B^Sz, and spins in a bistable anisotropy potential 
H. — —D S^ — B-S (supcrparamagncts) . The first one, with its cquispaccd spectrum, is 
a rotational counterpart of the quantum harmonic oscillator, while the anisotropic spin 
corresponds to problems of translational Brownian motion in non-harmonic potentials 
(double- well or periodic). The coupling to the bosonic bath considered have the 
structures Tish ~ r/ • (bilinear) and Tish ~ {Sz, (non- linear). The former 

may describe Kondo coupling to electron-hole excitations and the latter interaction 
with phonons, two important mechanisms in solids. Classically they correspond to 
field-type and anisotropy-type fluctuations in the spin Langevin equations. 

Both for isotropic and anisotropic spins we have given examples of static response, 
the dynamical susceptibility (to analyse the contribution of the different felaxation 
modes), and of spin resonance in transverse fields, which is very sensitive to the level 
spectrum and to the structure of the spin-bath coupling. Recall that effects like the 
spin resonance, or tunnel in transverse fields, demand the solution of the full density- 
matrix; such coherent dynamics involves off-diagonal elements and is not captured by 
a Pauli balance equation for the level populations. Finally, in some examples we used 
parameters close to actual qiiantum supcrparamagncts and typical experiments. 

We touched in passing the issue of the validity range of a master-equation 
description. Several limitations are inherited from the approximations required to 
derive quantum master equations (factorizing initial conditions. wc;ak system-bath 
coupling, high-T or semiclassical bath, etc.) along with manipulations specific of the 
problem addressed (secular or rotating-wave approximations, decoupling or adiabatic 
elimination of off-diagonal elements, etc.). This issue, however, is independent of the 
question of resolvability of quantum master equations by continued fractions, method 
which could in principle be applied to improved equations. 
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Appendix A. Energy-level spacings and effective fields 

In this appendix we calculate the level separations and related quantities (effective 
fields) for a uniaxial spin in a longitudinal field Tid = —DSl — BzSz (our unperturbed 
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Hamiltonian) . Then the eigenstates of are eigenstates of Hd too, i.e., Hdlrn) — 
£m|™)5 with e„i = —Dm^ — B^ra. The energy differences between these levels 

A„,„ = £„-e„i A„„i = -[£)(« + m) + S2] (n - m) , (A.l) 

appear in the unitary part of the Heisenberg equation for the Hubbard operators 
[Eq. 13)] and control the m ^ n transition rate Wn\m = W{Anm) [Eq. I|12|l ]. 

Level spacings. In the master equations considered only transitions between 
adjacent levels enter. For n — m±l the energy differences (jA.l|) give 

A, „,„.±i = ± [Di2m ± 1) + B,] . (A.2) 

In contrast with the equispaced Zeeman spectrum, the level spacings depend on m due 
to the anisotropy. For succesive pairs they are related by Am-i,m = Am,m+i — 2_D, 
decreasing as the barrier top is approached (Fig. QJ. To illustrate, for integer S at 
zero field the evolution from wells to barrier is 

Aw = As-i^s = D{2S - 1) ^ D{2S -3)^ > 3D ^ D ^ Ao,±i = At . (A.3) 

The boundaries coincide for 5 = 1, while Aw ^ 2DS for large S. For D ^ 0.5 K and 
S* = 10 (as in Mni2) we have the limit energy scales A„ ~ lOK and Ab ~ 0.5 K. 
Finally, when proceeding towards the classical limit fixing the anisotropy barrier D 
and Zeeman energy SB, the levels approach a continuum as A ~ 1/S. 

Effective, anisotropy, and critical fields. Classically one defines the effective field 
Bgg = —(dT-L/d^z), with /i^ = Scosd (the spin polar angle). This quantity enters 
in the Landau-Lifshitz precession equation. For a Hamiltonian including only the 
anisotropy term Ti, = —DS^, this definition gives the anisotropy field B^ — 2DS cosd. 
In the quantum case the /i^-derivative is naturally replaced by a finite difference 
Bcs = —(cm+i ~ £m)/^fn- Then, 6m = 1 plus Eq. I|A.2|) gives the effective and 
anisotropy fields of the quantum problem 

B^fi{m)=^D{2m + l) + Bz B^{m) = 2D{m + ^ . (A.4) 

For S ^ 1, we have i?a — 2DS{m/ S) 2DScosi3, recovering the classical result. 

To conclude, the critical field Be is that at which the barrier disappears (or 
equivalently the B^ that zeroes the last effective field i?eff(— 5); see Fig.^. Equating 
to zero the spacing between the last two levels, A^s+i.-s{Bc) = 0, one gets 

B, =D{2S-1), (-Balwelk). (A.5) 

This gives B^ — for S — 1/2, where there is no barrier, while the classical value 
matches the anisotropy field at the wells, 2DS = B^l^^o (~ 10 T in Mni2). 

Appendix B. Angular- momentum ladder factors 

Here we discuss some properties of the ladder factors = [S{S + 1) — m(m ± 1)]^^^. 
In the standard basis of eigenstates of and Sz, where S'^\m) = S{S + l)|m) 
and Sz\m) = m\m), the characterize the action of the raising/lowering operators 
S±\m) = if^\m ±1). In addition they are the expansion coefficients of S± ~ Sx ^ iSy 
on the Hubbard operator basis S± — J2m^m^m±i [EQ- 
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Figure Bl. Functions £{z) = Vl — z'-^ and £{z) = 2z\/l — z'^, sketching the 
dependence of the ladder factors Irn and Im on z ^ m/S. Right: squared factors 
£'^(z) = (1 — z^) and i'^{z) = {2z)-^{l — z^), as they enter in the relaxation term. 



Ordinary factors £rn- We introduce several alternative notations convenient in 
different contexts. The 2-index form is £m,m' — [S{S + 1) — mm']^^^, whence 

etr = ^™,™±i = VSiS+l)-m{m±l) = em±i,m ■ (B.l) 

This gives an explicit index connection in the ladder action S±\m) = im.,m±i\'m- i 1) 
as well as in S± = ^rn^rn,m±iX^^i. Both ira,m±i can be expressed in terms of a 
single ladder factor 

^m.m+l^im, ^m,m-l = ^m-l = S {S + I) - m{rn + 1) . (B.2) 

At m and — (to+ 1) the Im take equal values: Im — ^_„i-i- Thus is — ^-s-i — (end 
points), ^5-1 = t-S = V2S, £s-2 = £-s+i = ^2(25 - 1), . . . , £s-k = ^-(s+i)+fe = 
^^(25"+ 1 - k). Identifying m/S* ~ z = cost?, we have the behaviour ^ (1 — z^), 
which is the factor accounting for the reduction of the configuration space as the poles 
are approached fFig. IBl| . 

The "bar" factors £m- For the spin-bath coupling F ^ {Sz, S±} we come across 
some modulated ladder factors in the master equation: £m.m±i — (2m ± l)im,m±i- 
Using 2m ± 1 — m -\- (m ±1) we can write in a symmetric way 

£ra.m' — £m' .m — (m -\- m ^£m.m' ■ (B-3) 

Again we can introduce a compact single factor £m- 

£rn,m+l — £m 7 ^m,m — 1 £m~l : £m — {^rfl ^ . (B.4) 

The symmetry £^„i = £m-i, together with (— 2to + 1) = —(2m — 1), yields the 
corresponding (anti)symmetry of the bar factors, namely £_„i — —£m~i- Then 
is = -£~ s~i = ( boundaries), £s-i = -i-s = {2S - 1)^ 2^, 4-2 ^ -£ ~s+i = 
(23 - 3)^2(25-1), . . . , = e-{s+i)+k = (25 + 1 - 2k)^k{2S + 1 - k). 

For half- integer spin, £m vanishes at m = — 1/2. In general £m goes close to 
zero for small m (barrier top) as it has opposite signs for positive and negative m 
(Fig. IBll middle). For 5=1/2 all relevant £m vanish, reflecting that the coupling 
F ~ {Sz , S±} does not produce relaxation on a 5 = 1/2 spin; physically this F arises 
from the modulation of the anisotropy —DSz by the lattice vibrations 32 , but 5^ does 
not change the energy of 5 = 1/2 [recall the equilibrium result m^ = ^ th (^ Bz/T)]. 
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Appendix C. The Hubbard (level-shift) operators 

In this appendix we discuss some properties of the operators X™ = |n)(m| and derive 
their Heisenberg equations of motion in the conservative case. They form a complete 
set; if we think of a spin operator as a (25* + 1) x (25' + 1) matrix, then X™ is the 
matrix with zeroes everywhere, except one 1 at the position {n,m) [351 Ch. 1]. In this 
space of linear operators, they form an orthonormal system with respect to the scalar 
product (X, Y) EE Tt{X+Y). 

Properties of . We begin demonstrating several useful results. 

• Expressing an operator A in the Hubbard basis 

A = , A„™ = {n\A\m) . (C.l) 

Proof: using twice the closure relation I = ^^^^ \m){m\: 
A = i:»{n\AY.Jm){m\ = E^J'^I^I"^) I'^X"^! = E„™^n™ • 

• Equal-time relation 

X^xr = SkiX^ . (C.2) 
Proof: using the \m) basis orthonormality, X^Xf^ = |rt) (to| ~ 6kiX™. 

• Commutator 

[^n ; = ^klX"^ - ^nmXf . (C.3) 

Proof: [X^,Xp] = X^X^ - XpX^ plus the equal-time relation 

• Adjoint (recall the zeroes plus 1 at {n,m) representation of X^^) 

{x::^)^ = x-. (C.4) 

Proof: we use the auxiliary notation (-0, for the scalar product: 

The validity of this result V V & gives Eq. (|C.4p . 

• Relation with the density matrix (note the index ordering) 

(^r> = 9mn ■ (C.5) 

Proof: use the level-shift property X™\k) — 5mk\n) and (A) = Ti{gA): 
Tr (^.X™) = Efc(fc|e^"l^> = j:k{k\9\n)5mk = Qmn ■ 

As corollaries, replacing p — > I one gets the trace formula Tr (A™) = Snm- 
Then one can prove the orthonormality of the Hubbard basis: (A™, Xl) = 
Tr[(A™)+A[] = Tr(A;;A^) - Tr(<5„fcA^J - SnkSmi. 

Heisenberg equation of motion. As the Hubbard operators A™ do not depend 
explicitly on t, their Heisenberg dynamical equation is simply iA™ = [A™ ,7i]. Here 
we derive it explicitly for a Hamiltonian comprising a diagonal part 7id{Sz) (in the 
standard basis) plus a transverse field: H. = Ti.d{Sz) — ■^{B+S- -f B^S+). 

Let us commute the different parts of H with A™. For Hd we can write 
Hd = J2k^kX^, since its matrix elements are (n\Hd\m) = £„i5„m [see Eq. (IC.l|l ]. 
Then, we work [A™ ,Hd] out with help from Eq. (|tl3)l : 

[Xn , Tid] = J2k^k [A™ , A^'] = J2k^k{SmkX!^ - 5fe„A™) = -A„mA™ , 
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where A„,„ = e„ — Sm ■ For the transverse components we employ the representation 
S± = X]fe ^fe,fc±i^fe±i [Eq- O] and again commutator (|C.3|I : 

/; x^m=pl /J /^T V^^T^ 

where we have used im,m' = ^m',m and £m,m±i = ( [Appendix B| ). Finally gathering 
these partial results, including —B±/2, and multiplying across by — i 

i:r = -AK^ , Hd] + (i/2)i?+ [x;:' , + (i/2)i?_ [x;:' , 5+] , 

we arrive at the unitary evolution equation (|3J| for the Hubbard operators. 

If we had used the exact eigenstates of the total Hamiltonian, the term iA„mX™ 
would have been the only term in the equation (with B± entering via the £m)- Using 
the angular-momentum basis instead (for example, when the exact eigenstructure 
is not known or, for convenience, to treat time-dependent -B±), one needs to add 
explicitly the contribution of the transverse terms, as we have done here. 

Appendix D. The relaxation term 

In contact with the bath the closed equation of motion is augmented by the 
relaxation term ©. Here starting from such i?™ we derive the time-local relaxation 
term and then we particularise it to spin-bath couplings via S±. Finally we 
discuss the adoption of the secular or rotating-wave approximation and the issue of 
the positivity of the reduced description. 

Appendix D.l. Derivation of the time-local relaxation term Ml]) 

We start from ^R^,' = /!^dr{/C(r - i)F(r)[F,X™] - K{t - r) [F , X;"]F(r)}, 
where operators without argument are evaluated at t. The formal dependence on 
the previous evolution enters through F{t). We expand it on the Hubbard basis, 
— 'Yliki-FkiXl.{T), and replace the retarded dependences by the (dominant) 
5z-part of the conservative evolution X\^{t) — e~'*-*~'^-'^'''X^,(i), whence F{t) = 
Sfci -^fci^~'^*^^'''^'°'^A;(^)- Next, the change of variable s = t — t brings the integral 
relaxation term into the form 

-K' =E..^fc;^°°d.{/C(-s)e-'^^-Xi.[F,X:'] -/C(s)e"'^'^- [F,X:^]xI] . 

The s-dependences only occur in the kernel /C and the oscillating factors. On 
integrating them brings into scene the relaxation rates W^n = W{Aki), with W{A) = 
/o°°dse-^'*^/C(s) [Eq. lO]. Then, using /C(-t) = [/C(t)] * and Aik = -Am, gives 

- K = EHFki{w;^, xi[F,x:^] - . (d.i) 

This completes the task of getting a time-local relaxation term. Note however that 
i?™ still depends on the X in a non-explicit way (since F also contains them). % 

The rest of the calculation consists of simplifying the X structure of the above 
term. First, we expand the inner as F = J2k'i' Fk'i'Xl, and use [X^, = 

f It is assumed that at the initial time, to = — oo, system and bath were decoupled. If this is assumed 
to occur at some finite to, the rates adquire a time dependence W{A;t) = J^~^° ds e~^''^K{s). 
However, as in most problems system and bath had coexisted for a very long time, one shifts to — > — oo 
and the transients due to such decoupled initial conditions are washed away. 
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Sni'XJJ} — S„ik'Xl^ to work the commutators. Then we multiply by Xj. (from the left 
and right) and use the equal-time relation — (5;„X™ to get expressions linear 

in the Hubbards. Eventually, the result of the left and right multiplications is in turn 
multiplied by W^;*j,-Ffei and Wk\iFki respectively [Eq. (|D.l|l ] and summed over k i. 

Y.kl^l\kPkl Xl [F , X^] = J2kl^l\kPklFln Xl!" - J2kl'^n\kFknF„U' X^ 

Y.kiWk\iFki [F , Xi = Y.ik'W,n\iF„.iFk>n Xl - Y.kiWk\iFkiF..,k < • 

To clarify the structure and to arrive at the sought form X^n'm' ^™ri'" ^™ make 
several index changes. First line, first term: / I (summed index), k ^ n' , and we 
introduce X^m' ^mm'- First line, second term: /' — s- m' and k — > n' . Second line, first 
term: k' n' and I ^ m' . Second line, second term: k ^ i (summed index), I m' , 
and introduce (5„„' . Then we obtain for the above right-hand sides 

J2n'm' [^mm'{J2iW^^,Fn>iFen) - W*\^,Fn'nFmm']X^ 
^n'm' \y^m\m' Frriin' Fn' n ^nn' {^^£^^£[771' Firri' Fm^j'^ X.^/ . 

Substraction of these lines according to Eq. HU.1|I gives finally the relaxation term pif) . 
Appendix D.2. Relaxation term for couplings via S± 

Here we particularise the relaxation term to the coupling (|13|) . with F{S) 

"linear" in 5*^. We start expanding the anticommutators and introducing F± = 
ri^[v{Sz) S± -)- S±v{Sz)]: so that F ~ F^ + F^. To obtain the matrix elements 
Fnm — {n\F\m) we compute first those of using S-\m) = 4„.m-i|"i — 1)- 

(F„)„,„ = Lm,m-1 S„,m-1 , ^m.m' = ?7+[^^(™) + v{'m')]£m,rn' ■ (D-2) 

Note that the extended ladder factor Lm,m' is in general complex. Then by means of 
F+ = (F-)+ we arrive at = L„i,ra-i5n,7n-i + L*^+i^rn^n,m+i [i.e., Eq. lO]. 

Let us proceed to do the sums in the relaxation term with these matrix 
elements and some care. For the third line we need W^rn' FmeFim' X™ . The 

sum over m' gives 

Multiplying by F,„£ = Ljn+i,mSe^rn+i + ^m^m-i'^^.m-i and summing over i we obtain 

3rd line = L*j^,m-lLm~l,m-2^m-l\m~2 X'^ ^ 

+ ( \Lm,m-l\'^ W^m-l|m + |-t'm+l,m|^ W^m+l|m 
+ Lm+2,m+lLm+l,m W„i+i\m+2 ^r™^^ • 

The perturbative paths are clear, e.g., m + 2 ^ m-\-l mto connect X™"*"^ with X™. 
Inspection of Eq. reveals that the first line follows from the third by exchanging 
m ^ n, conjugating, and permuting upper and lower indices in X (adjointing): 

1st line = i„,„_iL„_i^„_2 W^,t-i|n-2 -'^"-2 
+ ( |L„,„_ipM^:_i|„ + |i„+i,„pM^:+i|„ )xr 

,7-* 7-* jrr* Y-m 

' ^n+2,n+l-^n+l,n n+l\n+2 ^n+2 ■ 

Now we are left with Y^n'm' i^n\n' +W^m|m')^n'n^™m' -'^t™' > t^c Central line of Eq. 
Exclude Fn'n and do first the sum over m', then multiply by Fn'n written in the form 
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Collecting the three contributions we finally obtain the specialisation of the Markovian 
relaxation term to the coupling (|13|l : 

+ |L„,„-ip W^„*_i|„ + ) X™ (D.3) 
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On invoking the secular approximation (see below) terms involving L x L ot L* x L* 
are dropped (last four lines) and only the terms of the type L x L* are retained (first 
four). This gives the secularized relaxation term (|15|l in the main text. 



Appendix D.3. The non-secular terms and the positivity issue 

Neglecting the "non-secular" terms L x L and L* x L* corresponds to the rotating- 
wave approximation of quantum optics, where counter-rotating terms are dropped. 
Specifically, writting schematically the coupling as Tish ~ S^{a~^ -\-ea )-\-S+{a -\-ea^), 
one can see that setting e = the last four lines of Eq. ljD.3|l dissappear To 
justify this manipulation, it is argued that 0+5*+ ^ e'^'^^'^)* oscillates faster than 
0+5- ~ gi('^-A)t g^j^j averaged out. A word of caution though. While this 

reasoning may apply to systems with monotonous spectrum (harmonic oscillator, 
isotropic spin), the justification is not satisfactory if the sign of Am,m,-i-i depends 
on m (upon A ^ — A it would be a'^ S- the term oscillating fast and a^5+, which is 
dropped, the slow one). 

Keeping all terms in Eq. (|D.3I) does not pose big problems for a continued-fraction 
solution. Only the line marked with the arrow (fifth one) breaks the 3-term recurrence 
in n, giving c„ ~ Q«,„-2C„-2 + Q«,n-ic„_i -|- Q„,„c„ + Q„,„+iC„+i + (Q)„,„+2C„+2, as 
equation of motion for (c„)„j = (X™). But this can be treated by 5 3 recurrence 
folding with block vectors and matrices 0] 1201 

Q — f '^2?i I Q ' ^ I Q2n+l,2n' \ 

" \ C2„+l J Y Q2n,2n' + 1 Q2n+l,2n' + l / ' 

recovering the canonical form C„ ^ Qn,n-iC„_i + Q„,„C„ + Q„,,i+iC„+i. 
Alternatively, one can handle c„+2 in an iterative way (keeping the ordinary vectors 
and matrices) which converges quickly for weak damping. Notwithstanding this, we 
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preferred to illustrate the continued-fraction treatment of the master equation with a 
simple yet generic case. 

Incidentally, the density-matrix equation without non-secular terms is of the so- 
called Lindblad type |^ for isotropic spins. In general, however, Eq. (jD.3|) does not 
ensure that for arbitrary initial conditions p is a positive operator at all times (which 
is required for its probabilistic interpretation). Although this sounds alarming, no real 
problem exists when one recalls the assumptions under which the master equation is 
derived and the meaning of g as partial trace of gtot- 

First, to derive the master equation one assumes that at some initial time to 
system and bath were decoupled £itot(^o) = Qifo) (g) Qhito)- The resulting evolution 
equation for g — Trb(f?tot) is precisely Eq. (|D.3p but with some time-dependent 
coefficients (recall the last footnote). But as this master equation comes from the 
solution of i(dptot/di) = [Ti. , gtot]) which preserves positivity of gtot (and hence of any 
partial trace), then g{t) > should hold to the accuracy considered in the system- 
bath coupling. Actually, no positivity violation has been reported when using the 
proper time-dependent coefficients. The cases oi g < explicitly shown |^ |^ 
correspond to use the master equation with the asymptotic i-independent coefficients 
(valid at long times), to address arbitrary initial condition problems. However, in the 
asymptotic regime, system and bath have been interacting for a long time and the 
possible g{t = 0) = Trb(£itot) cannot be chosen at our will. 

Except in problems like the atom inserted in an electromagnetic cavity or the like, 
decoupled initial conditions are just an mathematical trick to facilitate the obtaining 
of the form of the master equation. For problems like a spin in a solid, a Brownian 
particle in a fluid, etc., it is not natural to assume that system and bath just met, but 
that they were in contact for a long time. Thus, one shifts the initial time to — > —oo, so 
that at finite times (say t > Q) those factorizing initial conditions have been forgotten 
(except maybe in marginal cases) and system-plus-bath have evolved into some joint 
stationary state by their internal dynamics. Then the coefficients to be used are 
the asymptotic ones. But then it is natural that the initial conditions for g{t = 0) 
cannot be set arbitrarily, but only those compatible with the evolution of the coupled 
system-plus-bath. It is only overlooking this that the positivity becomes a issue. 

Finally, in the asymptotic t > regime, one can manipulate the system with fields 
or forces to prepare "initial conditions" for the problems of interest. This will only 
result in some time dependences of the master equation coefficients. In our case they 
arise when evolving X[s) back in time in the memory kernel, which now has to be 
done as X{s) = X{t) exp[— i f^d u A(m)], with the proper time-dependent Hamiltonian 
evolution ■ If the manipulation on the system is slow compared to the kernel width 
(bath correlation time) , one can use the leading term X{s) ~ X{t) e-i('*-*)A(t)_ Then 
one obtains precisely the ordinary master equation, with the asymptotic form of the 
coefficients, but parametrically modulated by the time-dependent fields or forces. 

Preparing the system in this way, the initial conditions follow from the master 
equation dynamics (and there is a single continuous process from to = — oo, instead 
of "initial conditions" being set twice). Then, to the accuracy in the coupling 
considered, the possitivity, hermiticity, normalization, etc., are preserved inasmuch 
as i(dgtot/dt) = \H , gtot] underlies the reduced dynamics for g = Trb(gtot)- 
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Appendix E. Transition and relaxation rates 

The rate function W^(A), associated to the kernel /C(t), can be expressed in terms of 
the bath spectral density [Eq. ((T?)|l ]. This follows integrating the relation (fTH)! between 
/C(t) and with Ke[J^dT exp(— irA)] — 7r^(A). In this appendix we consider 

W{A) for spectral densities of the form J{lj) — Xuj°'. We address two cases: (i) 
Ohmic bath a = 1 and (ii) the custom a = 3 super-Ohmic bath (photons or phonons 
in three dimensions). Note that in both cases a is odd. 



Unified functional form. Introducing the reduced variable u = A/T and omitting 
proportionality constants we write Eq. (fT?)|l as W{u > 0) ~ u°' Uu and W{u < 0) = 
+ 1) with the Bose factor Uu = l/(e" — 1). Let us first demonstrate that a 
single functional form gives the rates for both u < and w > 0: 



qu_i J i-e-u e-"-l 

Note that we have used the oddness of the spectral index a to absorb the sign change of 
the denominator (otherwise a sgn(u) appears). Thus W{—u) has the same functional 
dependence as the positive u case, i.e., (— u)" n^u, and we can write 

J(cj)cxtj" (odd a) W{u) = — Vu . (E.l) 

Figure IeTI displays the dependences oiW owu — A/T for a — \ and 3. For large and 
negative energy differences (emission of quanta) the behaviour is In the Ohmic 

case, the rate decreases monotonically from 1^(0) = 1 for positive energy differences 
(absorption). The super-Ohmic W{u), in contrast, after an initial parabolic take-off 
from zero, has a maximum at u 2\/2 and decreases exponentially dominated for 
large energy absorption. 

Derivatives ofW{u). These are needed in the perturbative calculations fSec. 17.21) . 
Having a unified functional form [Eq. ljE.l|) ] we can differentiate without worrying 
about distinguishing between absorption and emission processes: 

W {u) = + — . E.2 

^ ' e" - 1 (e" - l)(e-« - 1) ^ ' 

In the perturbative treatment of the probing field we actually need the i?2-derivatives 
of Wm\m±i- Applying the chain rule dW/dB^ — W'{u) du/dB^ we obtain a ± sign 
depending on the transition "lowering" or "raising" the second index [see Eq. (TOll ]: 

Wm\rn±l ^ A,n,m±l = ± [D{2m ± 1) + S^] dB,A.m,m±l = ±1 . 

Finally, one must recall to multiply Eq. l|E.2p also by 1/T since u — A/T. 



Relaxation rate. The combination T{u) — W{u) -f W{—u) appears is several 
problems (isotropic S — 1/2 and S — 1 spins, Sec. 18.3. H anisotropic S — 1, Sec. 19. 3.1)) . 
Using Eq. (|E.ip this symmetrized rate can be written as 

r = T4^(u)-hiy(-it) =u" cth(iu) . (E.3) 

This function behaves as 2u"~^ for small u, then grows monotonically and for large 
energy differences goes as u" (Fig. lEl| right). The monotony is proved differentiating 
Eq. l|fe3|l . whence F' = ^u°'-^{ashu - u)/ sh'^i^u) > 0, for a > 1. 
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Figure El. Transition rate W{u) (left panel) and relaxation rate F = 
W (u) -\-W {—u) (right). Curves for both the Ohmic bath, J(u)) oc lu (dashed lines), 
and the super-Ohmic bath J{uj) oc u)^ (solid) are shown. The limit functional 
dependences on u = A/T are marked. Introducing a high-frequency cuf-off ljy> 
in J(u)), both W(u) and r(u) would eventually drop to zero for \u\ 3> loyi/T, 
reflecting that no quanta of arbitrary large to are available to interchange energy 
with. 



Appendix F. Integral relaxation time 

The advantage of this quantity is that it can be obtained from the low frequency 
response x(^) — x(l ^ if^^int + ■ . .). This result follows from the definition Tint = 
J^At5M{t)/6M{Q) plus a low O expansion in the linear-response-theory relation 
X{n) =x(l-if^ j^dte-'^^[5M{t)/6M{id)]). Here we discuss a generalisation of the 
results for Tint of Refs. |21 OH by considering a generic discrete system obeying 
balance equations of the form |47| : 

Pm\m' is the m' ^ m transition probability, assumed to depend on the level spacings 
Pm\m' = P{^mm') and to obcy detailed balance Pm\m' = ^'^ Pm'\m- Besides 

the level separation is controlled by some bias F as A = A^^^ -I- F. 

Next one augments F by an oscillating probe F —f e'^*-f c.c. and seeks for a 

solution of Eq. (|KT|| of the form N,-n = + 5/ (A^m ^ e'"* -f c.c.) . The corresponding 

susceptibility is x(^) ^ X)™ (^)- 1*^^ frequencies the equations for Nm^ can 

be solved analytically and one gets ~^ 0)j whence Tint follows as |^ App. B] 

m ^m+l\m k=mi 

Here M' = dM/dy, with y F/T and Af = Em is the static response. The 

range of m is [mi, mt] and the superscript (0) indicates absence of probing field. Note 
that the auxiliary function <i>„i, in contrast to Pm|m'j depends only on equilibrium 
properties. 

To get Tint for isotropic spins, we simply cast the corresponding balance equations 
(Sec. I5.1|l into the form t^b\l\\ . by identifying Pm\m' = ^mm'^m\m'- To do so we 
have taken advantage of the 2-index form (jB.l|) of the ladder coefficients {£m-i = 
4n m-i = f-m-i m: and 4n = m+1 = m) , and recalled the simplyfied relations 
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Wm\m-i — a-nd VFm-i|m — Wo e ^ (Sec. |OJ. Thus Tint in this case is given by 

Tint = (1/mi) i'i'l/NjShlwo) , with = J:T=-s V. - fc) EH. 

To obtain Tjut for anisotropic spins we similarly bring their balance equations 
into the above form fSec. 19.3^1) . We use the 2-index notation of the modified factors 

= -^m.m+i and lm~i = ir,i,m-i plus lm,r,i' = jm' ,m [Eq. IIB.3|l] and identify 
Pm\m' = ^mm'^m\m'- Thcu the generic expression HF.2|) gives the relaxation time 

of the anisotropic-spin problem as Tint = I]m (*^m/-^™^^m^m|i|m)> with the 

same expression for (cf. Eq. (5.13) in Ref. and Eq. (16) in Ref. |3). RecaU 
finally that, in contrast to the isotropic case, no simplification occurs in the rates 
Wjn+i\m because of the non-equispaced spectrum. 



Appendix G. Continued-fraction methods to solve recurrence relations 

We conclude with a hands-on summary of the method of resolution of 3-term 
recurrences of the form 

Qia-i + Q.C. + Qfa+i =-f^, i = 1, 2, 3, . . . (G.l) 

Here the Qi are given coefficients {Q^ = 0) and the fi forcing or source terms. To 
obtain the unknown Ci one inserts in Eq. IjG.ip the ansatz 0] 

Ci — SiCi^i + Oi , (G.2) 

obtaining the following relations for the ladder coefficients Si and shifts a^: 

S.. = % , a, = - ^' + ^{;°'+^ . (G.3) 

Qi + Qi Si+i Qi + Qj Sij-i 

For finite recurrences Ciyi — for some /. To enforce this we set 5"/+! = 0, aj+i = 
and iterate downwards in ljG.3p getting all Si and down to i = 2. Now, to reconstruct 
all Ci from Eq. I|G.2|I . we only need the anchor Ci, which obeys: 

(Qi + 0^52) Ci = - (/i + Q+as) , (G.4) 

[Eq. (|G.1|) at i = 1 plus C2 — S2C1 +02]. Then, starting from the so-obtained Ci, we 
iterate Ci = SiCi-i -\- ai upwards, getting the solution of the recurrence (jG.ip . 

The above is the algorithmic form (easy to implement in a computer). The 

solution can also be written as Ci = (^Y[k=2 ^k^^i + J2k=2 (^TVj=k+i ^kj'^k- To 
illustrate the structure: C5 = S5S4S3S2C1 -\- S5S4S3a2 -\- S^S^a^ -f 5*5 04 -f 05. 
For homogeneous recurrences {fi = => = 0) the solution simply reads Ci = 
[SiSi-i ■ ■ ■ S2)Ci. 

As for the name of the method, note that Si is given in terms of a fraction with 
^i+i in the denominator, which can in turn be written as a fraction with 5^+2 in the 
denominator, and so on. This furnishes the structure of a continued fraction 

C^ ^^-p^- (G-5) 

qi H \ — 

92 H 

In simple problems one may identify the continued-fraction representation of some 
known function getting explicit analytical solutions. 

Differential recurrences like — Q~Ci-i -\- QiCi -\- Q^Ci^i -\- fi can be handled 
analogously for t-independent Qi. Laplace transformation plus g{s) = s g{s) — g{t — 0) 
brings the differential equation into the form Q^Ci-i + {Qi — s)Ci -t- Q^Ci+i = 



Continued- fraction methods for spin quantum-master equations 



41 



— [ft + Ci(0)]- Then, introducing some modified forcings fi = fi-\- Cj(0) and central 
coefficients Qi — Qi~s [cf. Eq. H36|) ] the equation acquires the structure of the ordinary 
recurrence relation (jG.l|) . and as such can be solved. § 

Further, the quantities in the recursions can be J- vectors (c^ and /J with 
J X J-matrix coefficients Qi. Then one proceeds analogously: inserting the ansatz 
Ci — SiCi-i + ai in the recurrence Q,^c,;_i + QiCi + Q^c^+i = —f^ and obtaining the 
coefficients and shifts a^, as in Eq. I|G.3|I . The only change is that the fraction bars 
stand now for matrix inversion ("from the left" A/B B^^ A), and one speaks of 
matrix continued fractions. The algorithm involves / inversions of J x J matrices, so 
reducing the storage requeriments and number of operations from the direct inversion 
(or diagonalisation) of the underlying (/ x J) x (/ x J) matrix problem. 

In the vector case, to iterate upwards Ci — §iCi_i + a^, the initial condition 
Ci obeys the matrix version of Eq. ljG.4|l . In the absence of forcing (/^ = ^ 
tti = 0; e.g., the 0th order Eq. l|23I)) the solution of such J x J system involves an 
overall multiplicative constant. One can add an extra equation to fix this, getting an 
augmented ( J+1) x J system |20l App. A] , which can be solved by methods appropriate 
for more equations than unknowns (e.g., Q^-decomposition) , yielding the required Ci. 
In our spin problems the normalisation condition 1 = ^^n = ^n(^n)n involves 
all c„ and cannot be used as the extra equation for c^^^s (which plays the role of 
Ci). However, we can fix arbitrarily one component, e.g., (c_s)_s — const (the extra 
equation) and normalise the solution at the end. A practical choice is a Gibbsian 
weight (c_5)_5 ~ Z"^exp(-e_s/T). 

To conclude, as the indices (n, m) can be half-integers (when S is so), we employ in 
the numerical implementation some integer indices i = n-\-{S-\-l) and j — m-l- (5*+ 1), 
running from 1 to 2S'-|-1 (= / = J). For integer S the equations can also be handled as 
two-sided recurrences — / < i < I (then the ansatz and initial conditions are slightly 
modified; see Refs. ^ 1201 App. A]). Although this may enhance stability is some 
problems, we have used the general protocol allowing for non-integer spins. 
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